#Code to generate final manuscript results
# Author: Ernani F. Choma
# Packages required: TAM, and, to plot Figure 2, sp, raster, and viridis.

# Please change the following lines:
# Line 8 -- working directory
# Lines 19-22, 27, 41, 46, 56, and 67 -- location of model inputs
# Lines 1027-1028 -- figure names and output location
# Line 2043 -- location of model input to generate supplementary dataset
# Line 2323 -- location of model input to generate figure 2

setwd("Results_with_Emissions/")

library("TAM")

###### 1. Reading Previous Outputs ######
#1.1. Marginal Damages -- outputs generated by the marginal damages model
#The marginal damages model generates outputs calculated with various CRFs, but we only use GEMM
load("Inputs/MV_Damages_BasePM_2017_Mortality_2017.RData")
load("Inputs/MV_Damages_BasePM_2017_Asthma_2019.RData")
load("Inputs/SRM_BasePM_2017_Mortality_2017.RData")
load("Inputs/SRM_BasePM_2017_Asthma_2019.RData")

#1.2. Emission Factors
#These are the emission factors for diesel school buses from EPAs NEI 2017, from Choma et al.
#link: https://www.pnas.org/doi/suppl/10.1073/pnas.2107402118/suppl_file/pnas.2107402118.sd03.csv
load("Inputs/School_Bus_Emissions_NEI.RData")


#These EFs include the TBW portion
#NEI 2017 used MOVES 2014, in which TBW emissions for Diesel School Buses are 0.0159 g PM2.5/mi, according to the EPA
#i.e., 0.0132 g/mi in brake wear + 0.0027 g/mi in tire wear
#Source: U.S. EPA https://nepis.epa.gov/Exe/ZyPDF.cgi/P100LCNE.PDF?Dockey=P100LCNE.PDF
#Creating separate variables for PM25 without TBW and for TBW emissions
DSB.EFs.NEI$PM25.Ex <- DSB.EFs.NEI$PM25_EF_grams_per_mile - 0.0159 #PM2.5 without TBW
DSB.EFs.NEI$PM25.TBW <- 0.0159 #TBW

#These are the GREET emission factors from Burnham (2021)
#Burnham, A. (2021). MOVES3 Vehicle Operation Emission Factors. 
#Available at: https://greet.es.anl.gov/files/update_moves3 (Accessed 7 August 2023)
load("Inputs/School_Bus_Emissions_GREET.RData")


#1.3. Grid mortality impacts, based on Choma et al. (2020, Environment International, https://doi.org/10.1016/j.envint.2020.106015)
#See a description of how the damages from Choma et al. (2020) were adjusted in the current PNAS paper -- see Methods section of main manuscript and supplementary information
load("Inputs/ESB_Deaths.RData")


#1.4. GHG benefits -- calculated using EPA's Social cost of Greenhouse Gases Application Workbook. For details, see main manuscript and SI
#Full citation for EPA's workbook
#U.S. Environmental Protection Agency (2024) 
#Supplementary Material for the Regulatory Impact Analysis for the Final Rulemaking, 
#“Standards of Performance for New, Reconstructed, and Modified Sources and Emissions Guidelines for Existing Sources: Oil and Natural Gas Sector Climate Review”
#EPA Report on the Social Cost of Greenhouse Gases: Estimates Incorporating Recent Scientific Advances: EPA Social Cost of Greenhouse Gases Application Workbook. 
#https://www.epa.gov/system/files/documents/2024-03/epa-sc-ghg-workbook_1.0.xlsx (Accessed 11 March 2024).
load("Inputs/GHG_BEN.RData")


#1.5. NCHS Urbanization classification: Source: CDC
# [dataset] U.S. Centers for Disease Control and Prevention, NCHS Urban-Rural Classification Scheme for Counties.
# Centers for Disease Control and Prevention: National Center for Health Statistics. 
# https://www.cdc.gov/nchs/data_access/urban_rural.htm (Accessed 15 June 2023).

# The exact file loaded on to NCHS.Urb can be downloaded at https://www.cdc.gov/nchs/data/data_acces_files/NCHSURCodes2013.xlsx
# Last downloaded on 04/20/24 at 1:45 PM
# It was opened in Excel and saved as .csv
NCHS.Urb <- read.csv("Inputs/NCHSURCodes2013.csv")
names(NCHS.Urb)
Urb <- subset(NCHS.Urb, select=c("FIPS.code","X2013.code"))
names(Urb) <- c("STCOU","NCHS.Urb")
str(Urb)
Urb <- Urb[Urb$STCOU %in% DSB.EFs.NEI$FIPS_STCOU,]
Urb <- Urb[order(Urb$STCOU,decreasing=FALSE),]

sum(Urb$STCOU == DSB.EFs.NEI$FIPS_STCOU)
sum(Urb$STCOU != DSB.EFs.NEI$FIPS_STCOU)


#1.6. School Bus Mileage parameters
#School Bus Mileage
#Source: Levinson et al. (2023), World Resources Institute
#M. Levinson, P. Burgoyne-Allen, A. Huntington, N. Hutchinson (2023) 
#Recommended total cost of ownership parameters for electric school buses: Summary of methods and data. Technical Note. 
#Washington, DC: World Resources Institute. 
#https://doi.org/10.46830/writn.22.00024 (Accessed 15 September 2023).
Year.VMT <- 14084
Total.VMT <- 14084*13.5

#Discounting future miles at 3% per year
aux.pv.mi <- (1/-0.03)*(exp(-0.03*13.5)-exp(0))
School.Bus.LT.Miles <- aux.pv.mi*(Total.VMT/13.5) 


###### 2. Adjusting Valuation ######
Asthma.VSC <- 610000 #This is 2022 USD (Appere et al., 2023)
# G. Appéré, D. Dussaux, A. Krupnick, M. Travers (2023) 
# Valuing a reduction in the risk and severity of asthma. 
# https://doi.org/10.1787/f289d29e-en (Accessed 20 June 2023).  

#2014 VSL
VSL.2014 <- 9.3 * 1e6
#Sources:
# U.S. Department of Health and Human Services (2016) 
#Guidelines for regulatory impact analysis. U.S. Department of Health and Human Services. 
#https://aspe.hhs.gov/sites/default/files/private/pdf/242926/HHS_RIAGuidance.pdf (accessed 7 November 2023).
# and
#L. A. Robinson, J. K. Hammitt, Valuing Reductions in Fatal Illness Risks:
#Implications of Recent Research. Health Economics 25, 1039-1052 (2016).


#Marginal Damages Model Output -- Adjusted to 2017 -- from Choma et al. (2021, PNAS, e2107402118, https://doi.org/10.1073/pnas.2107402118)
Income.Adjustment.2017 <- 351/334
Inflation.Adjustment.2017 <- 1.041166717 #World Bank 2015*2016*2017 GDP Deflator -- OLD FROM PNAS

#This study: Adjusting to 2022 USD and income levels
Income.Adjustment.2022 <- 362/334
Inflation.Adjustment.2022 <- 1.227406143 #World Bank 2015*2016*...*2021*2022 GDP Deflator


#Inflation adjustment from World Bank -- GDP Deflator
#https://data.worldbank.org/indicator/NY.GDP.DEFL.KD.ZG?locations=US 
#Accessed 23 Aug 2023

#Income Adjustment from BLS
#https://data.bls.gov/cgi-bin/srgate
#Series ID: LEU0252881600
#Accessed Aug 23, 2023

VSL.2017 <- VSL.2014 * Income.Adjustment.2017 * Inflation.Adjustment.2017 #Old Values from Marginal Damages Model
VSL.2022 <- VSL.2014 * Income.Adjustment.2022 * Inflation.Adjustment.2022 #New Values


#We also apply EPA's preferred cessation lag
#Source: 
# U.S. Environmental Protection Agency (2004) 
# Advisory Council on Clean Air Compliance analysis response to agency request on cessation lag. 
# https://nepis.epa.gov/Exe/ZyNET.exe/P100JMYX.TXT?ZyActionD=ZyDocument&Client=EPA&Index=2000+Thru+2005&Docs=&Query=&Time=&EndTime=&SearchMethod=1&TocRestrict=n&Toc=&TocEntry=&QField=&QFieldYear=&QFieldMonth=&QFieldDay=&IntQFieldOp=0&ExtQFieldOp=0&XmlQuery=&File=D%3A%5Czyfiles%5CIndex%20Data%5C00thru05%5CTxt%5C00000033%5CP100JMYX.txt&User=ANONYMOUS&Password=anonymous&SortMethod=h%7C-&MaximumDocuments=1&FuzzyDegree=0&ImageQuality=r75g8/r75g8/x150y150g16/i425&Display=hpfr&DefSeekPage=x&SearchBack=ZyActionL&Back=ZyActionS&BackDesc=Results%20page&MaximumPages=1&ZyEntry=1&SeekPage=x&ZyPURL# 
# (Accessed 4 August 2023).

Cessation.Lag <- c(0.3,rep(0.5/4,4),rep(0.2/15,15))
Discount.Rate <- 0.03
Year <- 1:20
Present.Value <- (1/Discount.Rate)*(exp(-Discount.Rate*(Year-1))-exp(-Discount.Rate*(Year)))
Sum.Present.Value <- t(Cessation.Lag) %*% Present.Value
Present.Value.Multiplier <- as.numeric(Sum.Present.Value)

Present.Value.Multiplier #Our marginal damages are multiplied by ~0.89.
#Undiscounted values can be obtained by removing this step, i.e. setting the Present.Value.Multiplier to 1
#Present.Value.Multiplier <- 1 #For undiscounted damages

VSL.2017.Adj <- Present.Value.Multiplier * VSL.2017
VSL.2022.Adj <- Present.Value.Multiplier * VSL.2022

# VSL.Adj.Factor <- VSL.2022/VSL.2017 #To adjust Marginal Damages Model outputs
VSL.Adj.Factor <- VSL.2022.Adj/VSL.2017.Adj #To adjust Marginal Damages Model outputs

#Adjusting Marginal Damages Model Results
Marg.Damages.GEMM_BasePM_2017_Mortality_2017 <- Marg.Damages.GEMM_BasePM_2017_Mortality_2017 * VSL.Adj.Factor
SR.Damages.GEMM_BasePM_2017_Mortality_2017 <- SR.Damages.GEMM_BasePM_2017_Mortality_2017 * VSL.Adj.Factor


###### 3. Spatial variation in emission factors ######
aux.probs <- c(0.025,0.05,0.25,0.5,0.75,0.95,0.975)

#Emission Quantiles
Emis.Distrib <- as.data.frame(rep("a",12))
names(Emis.Distrib) <- "EFs"
aux.polsem <- c("PM25","SO2","NOX","NH3","VOC","PMEX")
Emis.Distrib$EFs <- c(paste0(aux.polsem,"_Unweighted"),paste0(aux.polsem,"_Weighted_by_VMT"))
Emis.Distrib$pctl_025 <- 0
Emis.Distrib$pctl_05 <- 0
Emis.Distrib$pctl_25 <- 0
Emis.Distrib$pctl_50 <- 0
Emis.Distrib$pctl_75 <- 0
Emis.Distrib$pctl_95 <- 0
Emis.Distrib$pctl_975 <- 0
Emis.Distrib$mean <- 0
Emis.Distrib$sd <- 0

names(DSB.EFs.NEI)
aux.pols <- c(7:11,12)

for (i in 1:6)
{
  aux.col <- aux.pols[i]
  #Unweighted
  Emis.Distrib[i,2:8] <- as.numeric(quantile(DSB.EFs.NEI[,aux.col],probs=aux.probs))
  Emis.Distrib[i,9] <- mean(DSB.EFs.NEI[,aux.col])
  Emis.Distrib[i,10] <- sd(DSB.EFs.NEI[,aux.col])
  #Weighted by VMT
  Emis.Distrib[(i+6),2:8] <- as.numeric(weighted_quantile(DSB.EFs.NEI[,aux.col],w=DSB.EFs.NEI$VMT,probs=aux.probs))
  Emis.Distrib[(i+6),9] <- as.numeric(weighted_mean(DSB.EFs.NEI[,aux.col],w=DSB.EFs.NEI$VMT))
  Emis.Distrib[(i+6),10] <- as.numeric(weighted_sd(DSB.EFs.NEI[,aux.col],w=DSB.EFs.NEI$VMT))
}

#First row of Table 1
cbind(Emis.Distrib[c(12,8:11),1],round(Emis.Distrib[c(12,8:11),c(9,2,8)],digits=3))

#Rows 2-4 of Table 1
DSB.EFs.GREET

###### 4. Calculating impacts ######
#Creating data frame
DSB.Mi <- as.data.frame(matrix(rep(0,3108*36*5),ncol=36*5,byrow=TRUE))
#Rows are 1 for each county
#Columns represent combination of Emission factor ; Pollutant; & Outcome

#Outcomes are: 
#In USD: Deaths, Asthma, Both (Deaths+Asthma) 
#Attributable number: AttribD (attributable deaths) ; AttribA (attributable new childhood asthma cases)
Outcome.Names <- c("Deaths","Asthma","BothOutcomes","AttribD","AttribA")

#Pollutants
Pollutant.Names <- c("PM25","SO2","NOX","NH3","VOC","PMDTOX","PMTBW","SUM","SUMDTOX")

#Sets of Emission factors used
EF.Names <- c("NEI","MY2005","MY2010","MY2020")

aux.names <- expand.grid(Pollutant.Names,EF.Names,Outcome.Names)
names(DSB.Mi) <- paste(aux.names[,1],aux.names[,2],aux.names[,3],sep="_")

###### 4.1. Calculating per mile impacts of DSBs ######
#Mortality 
#NEI
DSB.Mi$PM25_NEI_Deaths <- DSB.EFs.NEI$PM25.Ex * 1e-6 * Marg.Damages.GEMM_BasePM_2017_Mortality_2017[,1]
DSB.Mi$SO2_NEI_Deaths <- DSB.EFs.NEI$SO2_EF_grams_per_mile * 1e-6 * Marg.Damages.GEMM_BasePM_2017_Mortality_2017[,2]
DSB.Mi$NOX_NEI_Deaths <- DSB.EFs.NEI$NOX_EF_grams_per_mile * 1e-6 * Marg.Damages.GEMM_BasePM_2017_Mortality_2017[,3]
DSB.Mi$NH3_NEI_Deaths <- DSB.EFs.NEI$NH3_EF_grams_per_mile * 1e-6 * Marg.Damages.GEMM_BasePM_2017_Mortality_2017[,4]
DSB.Mi$VOC_NEI_Deaths <- DSB.EFs.NEI$VOC_EF_grams_per_mile * 1e-6 * Marg.Damages.GEMM_BasePM_2017_Mortality_2017[,5]
DSB.Mi$SUM_NEI_Deaths <- DSB.Mi$PM25_NEI_Deaths + DSB.Mi$SO2_NEI_Deaths + DSB.Mi$NOX_NEI_Deaths + DSB.Mi$NH3_NEI_Deaths + DSB.Mi$VOC_NEI_Deaths
DSB.Mi$PMDTOX_NEI_Deaths <- 5*DSB.Mi$PM25_NEI_Deaths
DSB.Mi$SUMDTOX_NEI_Deaths <- DSB.Mi$PMDTOX_NEI_Deaths + DSB.Mi$SO2_NEI_Deaths + DSB.Mi$NOX_NEI_Deaths + DSB.Mi$NH3_NEI_Deaths + DSB.Mi$VOC_NEI_Deaths
DSB.Mi$PMTBW_NEI_Deaths <- DSB.EFs.NEI$PM25.TBW * 1e-6 * Marg.Damages.GEMM_BasePM_2017_Mortality_2017[,1]

#MY 2005
DSB.Mi$PM25_MY2005_Deaths <- DSB.EFs.GREET$PM25.Ex[1] * 1e-6 * Marg.Damages.GEMM_BasePM_2017_Mortality_2017[,1]
DSB.Mi$NOX_MY2005_Deaths <- DSB.EFs.GREET$NOX[1] * 1e-6 * Marg.Damages.GEMM_BasePM_2017_Mortality_2017[,3]
DSB.Mi$VOC_MY2005_Deaths <- DSB.EFs.GREET$VOC[1] * 1e-6 * Marg.Damages.GEMM_BasePM_2017_Mortality_2017[,5]
DSB.Mi$SUM_MY2005_Deaths <- DSB.Mi$PM25_MY2005_Deaths + DSB.Mi$NOX_MY2005_Deaths + DSB.Mi$VOC_MY2005_Deaths
DSB.Mi$PMDTOX_MY2005_Deaths <- 5*DSB.Mi$PM25_MY2005_Deaths
DSB.Mi$SUMDTOX_MY2005_Deaths <- DSB.Mi$PMDTOX_MY2005_Deaths + DSB.Mi$NOX_MY2005_Deaths + DSB.Mi$VOC_MY2005_Deaths
DSB.Mi$PMTBW_MY2005_Deaths <- DSB.EFs.GREET.TBW$PM25.TBW[1] * 1e-6 * Marg.Damages.GEMM_BasePM_2017_Mortality_2017[,1]

#MY 2010
DSB.Mi$PM25_MY2010_Deaths <- DSB.EFs.GREET$PM25.Ex[2] * 1e-6 * Marg.Damages.GEMM_BasePM_2017_Mortality_2017[,1]
DSB.Mi$NOX_MY2010_Deaths <- DSB.EFs.GREET$NOX[2] * 1e-6 * Marg.Damages.GEMM_BasePM_2017_Mortality_2017[,3]
DSB.Mi$VOC_MY2010_Deaths <- DSB.EFs.GREET$VOC[2] * 1e-6 * Marg.Damages.GEMM_BasePM_2017_Mortality_2017[,5]
DSB.Mi$SUM_MY2010_Deaths <- DSB.Mi$PM25_MY2010_Deaths + DSB.Mi$NOX_MY2010_Deaths + DSB.Mi$VOC_MY2010_Deaths
DSB.Mi$PMDTOX_MY2010_Deaths <- 5*DSB.Mi$PM25_MY2010_Deaths
DSB.Mi$SUMDTOX_MY2010_Deaths <- DSB.Mi$PMDTOX_MY2010_Deaths + DSB.Mi$NOX_MY2010_Deaths + DSB.Mi$VOC_MY2010_Deaths
DSB.Mi$PMTBW_MY2005_Deaths <- DSB.EFs.GREET.TBW$PM25.TBW[2] * 1e-6 * Marg.Damages.GEMM_BasePM_2017_Mortality_2017[,1]

#MY2020
DSB.Mi$PM25_MY2020_Deaths <- DSB.EFs.GREET$PM25.Ex[3] * 1e-6 * Marg.Damages.GEMM_BasePM_2017_Mortality_2017[,1]
DSB.Mi$NOX_MY2020_Deaths <- DSB.EFs.GREET$NOX[3] * 1e-6 * Marg.Damages.GEMM_BasePM_2017_Mortality_2017[,3]
DSB.Mi$VOC_MY2020_Deaths <- DSB.EFs.GREET$VOC[3] * 1e-6 * Marg.Damages.GEMM_BasePM_2017_Mortality_2017[,5]
DSB.Mi$SUM_MY2020_Deaths <- DSB.Mi$PM25_MY2020_Deaths + DSB.Mi$NOX_MY2020_Deaths + DSB.Mi$VOC_MY2020_Deaths
DSB.Mi$PMDTOX_MY2020_Deaths <- 5*DSB.Mi$PM25_MY2020_Deaths
DSB.Mi$SUMDTOX_MY2020_Deaths <- DSB.Mi$PMDTOX_MY2020_Deaths + DSB.Mi$NOX_MY2020_Deaths + DSB.Mi$VOC_MY2020_Deaths
DSB.Mi$PMTBW_MY2005_Deaths <- DSB.EFs.GREET.TBW$PM25.TBW[3] * 1e-6 * Marg.Damages.GEMM_BasePM_2017_Mortality_2017[,1]

#Asthma
DSB.Mi$PM25_NEI_Asthma <- DSB.EFs.NEI$PM25.Ex * 1e-6 * Marg.Damages.Tetreault_BasePM_2017_Asthma_2019[,16]
DSB.Mi$SO2_NEI_Asthma <- DSB.EFs.NEI$SO2_EF_grams_per_mile * 1e-6 * Marg.Damages.Tetreault_BasePM_2017_Asthma_2019[,17]
DSB.Mi$NOX_NEI_Asthma <- DSB.EFs.NEI$NOX_EF_grams_per_mile * 1e-6 * Marg.Damages.Tetreault_BasePM_2017_Asthma_2019[,18]
DSB.Mi$NH3_NEI_Asthma <- DSB.EFs.NEI$NH3_EF_grams_per_mile * 1e-6 * Marg.Damages.Tetreault_BasePM_2017_Asthma_2019[,19]
DSB.Mi$VOC_NEI_Asthma <- DSB.EFs.NEI$VOC_EF_grams_per_mile * 1e-6 * Marg.Damages.Tetreault_BasePM_2017_Asthma_2019[,20]
DSB.Mi$SUM_NEI_Asthma <- DSB.Mi$PM25_NEI_Asthma + DSB.Mi$SO2_NEI_Asthma + DSB.Mi$NOX_NEI_Asthma + DSB.Mi$NH3_NEI_Asthma + DSB.Mi$VOC_NEI_Asthma
DSB.Mi$PMDTOX_NEI_Asthma <- 5*DSB.Mi$PM25_NEI_Asthma
DSB.Mi$SUMDTOX_NEI_Asthma <- DSB.Mi$PMDTOX_NEI_Asthma + DSB.Mi$SO2_NEI_Asthma + DSB.Mi$NOX_NEI_Asthma + DSB.Mi$NH3_NEI_Asthma + DSB.Mi$VOC_NEI_Asthma
DSB.Mi$PMTBW_NEI_Asthma <- DSB.EFs.NEI$PM25.TBW * 1e-6 * Marg.Damages.Tetreault_BasePM_2017_Asthma_2019[,16]


#MY2005
DSB.Mi$PM25_MY2005_Asthma <- DSB.EFs.GREET$PM25.Ex[1] * 1e-6 * Marg.Damages.Tetreault_BasePM_2017_Asthma_2019[,16]
DSB.Mi$NOX_MY2005_Asthma <- DSB.EFs.GREET$NOX[1] * 1e-6 * Marg.Damages.Tetreault_BasePM_2017_Asthma_2019[,18]
DSB.Mi$VOC_MY2005_Asthma <- DSB.EFs.GREET$VOC[1] * 1e-6 * Marg.Damages.Tetreault_BasePM_2017_Asthma_2019[,20]
DSB.Mi$SUM_MY2005_Asthma <- DSB.Mi$PM25_MY2005_Asthma + DSB.Mi$NOX_MY2005_Asthma + DSB.Mi$VOC_MY2005_Asthma
DSB.Mi$PMDTOX_MY2005_Asthma <- 5*DSB.Mi$PM25_MY2005_Asthma
DSB.Mi$SUMDTOX_MY2005_Asthma <- DSB.Mi$PMDTOX_MY2005_Asthma + DSB.Mi$NOX_MY2005_Asthma + DSB.Mi$VOC_MY2005_Asthma
DSB.Mi$PMTBW_MY2005_Asthma <- DSB.EFs.GREET.TBW$PM25.TBW[1] * 1e-6 * Marg.Damages.Tetreault_BasePM_2017_Asthma_2019[,16]


#MY2010
DSB.Mi$PM25_MY2010_Asthma <- DSB.EFs.GREET$PM25.Ex[2] * 1e-6 * Marg.Damages.Tetreault_BasePM_2017_Asthma_2019[,16]
DSB.Mi$NOX_MY2010_Asthma <- DSB.EFs.GREET$NOX[2] * 1e-6 * Marg.Damages.Tetreault_BasePM_2017_Asthma_2019[,18]
DSB.Mi$VOC_MY2010_Asthma <- DSB.EFs.GREET$VOC[2] * 1e-6 * Marg.Damages.Tetreault_BasePM_2017_Asthma_2019[,20]
DSB.Mi$SUM_MY2010_Asthma <- DSB.Mi$PM25_MY2010_Asthma + DSB.Mi$NOX_MY2010_Asthma + DSB.Mi$VOC_MY2010_Asthma
DSB.Mi$PMDTOX_MY2010_Asthma <- 5*DSB.Mi$PM25_MY2010_Asthma
DSB.Mi$SUMDTOX_MY2010_Asthma <- DSB.Mi$PMDTOX_MY2010_Asthma + DSB.Mi$NOX_MY2010_Asthma + DSB.Mi$VOC_MY2010_Asthma
DSB.Mi$PMTBW_MY2010_Asthma <- DSB.EFs.GREET.TBW$PM25.TBW[2] * 1e-6 * Marg.Damages.Tetreault_BasePM_2017_Asthma_2019[,16]

#MY2020
DSB.Mi$PM25_MY2020_Asthma <- DSB.EFs.GREET$PM25.Ex[3] * 1e-6 * Marg.Damages.Tetreault_BasePM_2017_Asthma_2019[,16]
DSB.Mi$NOX_MY2020_Asthma <- DSB.EFs.GREET$NOX[3] * 1e-6 * Marg.Damages.Tetreault_BasePM_2017_Asthma_2019[,18]
DSB.Mi$VOC_MY2020_Asthma <- DSB.EFs.GREET$VOC[3] * 1e-6 * Marg.Damages.Tetreault_BasePM_2017_Asthma_2019[,20]
DSB.Mi$SUM_MY2020_Asthma <- DSB.Mi$PM25_MY2020_Asthma + DSB.Mi$NOX_MY2020_Asthma + DSB.Mi$VOC_MY2020_Asthma
DSB.Mi$PMDTOX_MY2020_Asthma <- 5*DSB.Mi$PM25_MY2020_Asthma
DSB.Mi$SUMDTOX_MY2020_Asthma <- DSB.Mi$PMDTOX_MY2020_Asthma + DSB.Mi$NOX_MY2020_Asthma + DSB.Mi$VOC_MY2020_Asthma
DSB.Mi$PMTBW_MY2020_Asthma <- DSB.EFs.GREET.TBW$PM25.TBW[3] * 1e-6 * Marg.Damages.Tetreault_BasePM_2017_Asthma_2019[,16]

names(DSB.Mi)
DSB.Mi[,(36*2) + 1:36] <- DSB.Mi[,1:36] + DSB.Mi[,36+1:36]
DSB.Mi[,(36*3) + 1:36] <- DSB.Mi[,1:36] * (1/VSL.2022.Adj)
DSB.Mi[,(36*4) + 1:36] <- DSB.Mi[,36+1:36] * (1/Asthma.VSC)


###### 4.2. Calculating per bus impacts of DSBs ######
grep("Attrib",names(DSB.Mi))
DSB.BUS <- DSB.Mi * School.Bus.LT.Miles
DSB.BUS[,c(((36*3) + 1:36),((36*4) + 1:36))] <- DSB.Mi[,c(((36*3) + 1:36),((36*4) + 1:36))] * Total.VMT

sum(c(((36*3) + 1:36),((36*4) + 1:36)) == 109:180)
sum(c(((36*3) + 1:36),((36*4) + 1:36)) != 109:180)


###### 4.3. Aggregating by urbanization
#Creating data frame
DSB.Mi.URB <- as.data.frame(matrix(rep(0,8*36*5),ncol=36*5,byrow=TRUE))
names(DSB.Mi.URB) <- names(DSB.Mi)

for (i in 1:7)
{
  Urb.cou <- which(Urb$NCHS.Urb == i)
  if (i==7)
  {
    Urb.cou <- which(Urb$NCHS.Urb %in% c(1,2))
  }
  aux.DSB.Mi <- DSB.Mi[Urb.cou,]
  aux.VMT <- DSB.EFs.NEI$VMT[Urb.cou]
  for (j in 1:dim(aux.DSB.Mi)[2])
  {
    DSB.Mi.URB[i,j] <- sum(aux.DSB.Mi[,j]*aux.VMT)/sum(aux.VMT)
  }
  
}

for (j in 1:dim(DSB.Mi)[2])
{
  DSB.Mi.URB[8,j] <- sum(DSB.Mi[,j]*DSB.EFs.NEI$VMT)/sum(DSB.EFs.NEI$VMT)
}


aux.which.attrib <- grep("Attrib",names(DSB.Mi.URB))
aux.which.attrib <- aux.which.attrib[order(aux.which.attrib, decreasing=FALSE)]
DSB.BUS.URB <- DSB.Mi.URB * School.Bus.LT.Miles
DSB.BUS.URB[,aux.which.attrib] <- DSB.Mi.URB[,aux.which.attrib] * Total.VMT



###### 4.4. Calculating ESB impacts
ESB.Mi <- as.data.frame(matrix(rep(-1,15),ncol=15))

#Pollutants
ESB.Pollutant.Names <- c("SO2","NOX","SUM")

aux.names <- expand.grid(ESB.Pollutant.Names,Outcome.Names)
names(ESB.Mi) <- paste(aux.names[,1],aux.names[,2],sep="_")
#Columns follow the same names as DSB.BUS -- except that only two pollutants are considered for ESBs/Power plants (NOX and SO2)
#We also assume that TBW impacts are the same for DSBs and ESBs

#Deaths
ESB.Mi$SO2_Deaths <- ESB.Deaths$USD.Mi[1]
ESB.Mi$NOX_Deaths <- ESB.Deaths$USD.Mi[2]
ESB.Mi$SUM_Deaths <- ESB.Deaths$USD.Mi[3]

#Accounting for Asthma damages -- assuming ratio deaths/asthma damages is the same for DSBs and ESBs
aux.Both.Avg <- sum(DSB.BUS$SUM_NEI_BothOutcomes * DSB.EFs.NEI$VMT)/sum(DSB.EFs.NEI$VMT)
aux.Deaths.Avg <- sum(DSB.BUS$SUM_NEI_Deaths * DSB.EFs.NEI$VMT)/sum(DSB.EFs.NEI$VMT)

Adj.Asthma <- aux.Both.Avg/aux.Deaths.Avg

#Sum
ESB.Mi$SO2_BothOutcomes <- ESB.Mi$SO2_Deaths * Adj.Asthma
ESB.Mi$NOX_BothOutcomes <- ESB.Mi$NOX_Deaths * Adj.Asthma
ESB.Mi$SUM_BothOutcomes <- ESB.Mi$SUM_Deaths * Adj.Asthma

#Asthma
ESB.Mi$SO2_Asthma <- ESB.Mi$SO2_BothOutcomes - ESB.Mi$SO2_Deaths
ESB.Mi$NOX_Asthma <- ESB.Mi$NOX_BothOutcomes - ESB.Mi$NOX_Deaths
ESB.Mi$SUM_Asthma <- ESB.Mi$SUM_BothOutcomes - ESB.Mi$SUM_Deaths

#Attributable deaths
ESB.Mi$SO2_AttribD <- ESB.Mi$SO2_Deaths * (1/VSL.2022.Adj)
ESB.Mi$NOX_AttribD <- ESB.Mi$NOX_Deaths * (1/VSL.2022.Adj)
ESB.Mi$SUM_AttribD <- ESB.Mi$SUM_Deaths * (1/VSL.2022.Adj)

#Attributable new childhood asthma cases
ESB.Mi$SO2_AttribA <- ESB.Mi$SO2_Asthma * (1/Asthma.VSC)
ESB.Mi$NOX_AttribA <- ESB.Mi$NOX_Asthma * (1/Asthma.VSC)
ESB.Mi$SUM_AttribA <- ESB.Mi$SUM_Asthma * (1/Asthma.VSC)

grep("Attrib",names(ESB.Mi))
aux.which.attrib <- grep("Attrib",names(ESB.Mi))
aux.which.attrib <- aux.which.attrib[order(aux.which.attrib,decreasing=FALSE)]

ESB.BUS <- ESB.Mi * School.Bus.LT.Miles
ESB.BUS[,aux.which.attrib] <- ESB.Mi[,aux.which.attrib] * Total.VMT

###### 4.5 Creating datasets presented in Figures and SI Figures ######
#Urbanization categories
Urb.Categories <- paste0("NCHS",c(1:6,"US"))

#Climate Benefits:
#a) per bus in Figs 1, S10, S11
#b) per mile in Fig 3 & 4
BEN.BUS.GHG <- GHG.BEN$ESB.BEN.BUS
BEN.Mi.GHG <- GHG.BEN$ESB.BEN.Mi

GHG.BEN$ESB.BEN.BUS/GHG.BEN$ESB.BEN.Mi

###### 4.5.1. Results for Figure 1 ######
#a) NEI EFs: Benefits per bus, by outcome & urbanization
#b) NEI EFs: Health impacts of DSBs per bus, by pollutant & urbanization 
#c) Health impacts of ESBs per bus, by pollutant
Fig1.BEN.BUS <- subset(DSB.BUS.URB, select=c("SUM_NEI_Deaths","SUM_NEI_Asthma")) #DSB Impacts
Fig1.BEN.BUS$SUM_NEI_Deaths <- Fig1.BEN.BUS$SUM_NEI_Deaths - ESB.BUS$SUM_Deaths
Fig1.BEN.BUS$SUM_NEI_Asthma <- Fig1.BEN.BUS$SUM_NEI_Asthma - ESB.BUS$SUM_Asthma
Fig1.BEN.BUS <- Fig1.BEN.BUS[-7,]

aux.rownames <- paste("SUM","NEI",Urb.Categories,sep="_")

Fig1.BEN.BUS <- cbind(aux.rownames,Fig1.BEN.BUS)
colnames(Fig1.BEN.BUS) <- c("Data_Name","Deaths","Asthma")


Fig1.DSB.BUS <- subset(DSB.BUS.URB, select=c("PM25_NEI_BothOutcomes",
                                             "SO2_NEI_BothOutcomes",
                                             "NOX_NEI_BothOutcomes",
                                             "NH3_NEI_BothOutcomes",
                                             "VOC_NEI_BothOutcomes"))
Fig1.DSB.BUS <- Fig1.DSB.BUS[-7,]
aux.rownames <- paste("NEI","BothOutcomes",Urb.Categories,sep="_")
Fig1.DSB.BUS <- cbind(aux.rownames,Fig1.DSB.BUS)
colnames(Fig1.DSB.BUS) <- c("Data_Name","PM25","SO2","NOX","NH3","VOC")


Fig1.ESB.BUS <- subset(ESB.BUS, select=c("SO2_BothOutcomes","NOX_BothOutcomes"))
aux.rownames <- paste0("BothOutcomes")
Fig1.ESB.BUS <- cbind(aux.rownames,Fig1.ESB.BUS)
colnames(Fig1.ESB.BUS) <- c("Data_Name","SO2","NOX")

#Showing values
cbind(Fig1.BEN.BUS[,1],round(rowSums(Fig1.BEN.BUS[,-1]),digits=-2))
cbind(Fig1.DSB.BUS[,1],round(rowSums(Fig1.DSB.BUS[,-1]),digits=-2))
cbind(Fig1.ESB.BUS[,1],round(rowSums(Fig1.ESB.BUS[,-1]),digits=-2))



###### 4.5.2. Results for Figure 2 ######
#NEI EFs: Benefits per bus, by county
Fig2.BEN.BUS <- DSB.BUS$SUM_NEI_BothOutcomes - ESB.BUS$SUM_BothOutcomes

###### 4.5.3. Results for Figure 3 ######
#a) ALL EFs: Benefits per mile, by outcome & urbanization
#b) ALL EFs: Health impacts of DSBs per mile, by pollutant & urbanization 
#c) ESBs: Health impacts of ESBs per mile, by pollutant

aux.ef <- c("NEI","MY2005","MY2010","MY2020")
aux.oc <- c("Deaths","Asthma")

for (i in 1:length(aux.ef))
{
  aux.subset <- paste("SUM",aux.ef[i],aux.oc,sep="_")
  aux.subset2 <- subset(DSB.Mi.URB, select=aux.subset)
  aux.subset2 <- aux.subset2[-7,]
  aux.d <- grep("Deaths",colnames(aux.subset2))
  aux.subset2[,aux.d] <- aux.subset2[,aux.d] - ESB.Mi$SUM_Deaths
  aux.a <- grep("Asthma",colnames(aux.subset2))
  aux.subset2[,aux.a] <- aux.subset2[,aux.a] - ESB.Mi$SUM_Asthma
  aux.names <- paste("SUM",aux.ef[i],Urb.Categories,sep="_") 
  aux.subset2 <- cbind(aux.names,aux.subset2)
  colnames(aux.subset2) <- c("Data_Name",aux.oc)
  if (i==1)
  {
    Fig3.BEN.Mi <- aux.subset2
  } else {
    Fig3.BEN.Mi <- rbind(Fig3.BEN.Mi,aux.subset2)
  }
}

aux.ef <- c("NEI","MY2005","MY2010","MY2020")
aux.pol <- c("PM25","SO2","NOX","NH3","VOC")

for (i in 1:length(aux.ef))
{
  
  aux.subset <- paste(aux.pol,aux.ef[i],"BothOutcomes",sep="_")
  aux.subset2 <- subset(DSB.Mi.URB, select=aux.subset)
  aux.subset2 <- aux.subset2[-7,]
  aux.names <- paste(aux.ef[i],"BothOutcomes",Urb.Categories,sep="_") 
  aux.subset2 <- cbind(aux.names,aux.subset2)
  colnames(aux.subset2) <- c("Data_Name",aux.pol)
  if (i==1)
  {
    Fig3.DSB.Mi <- aux.subset2
  } else {
    Fig3.DSB.Mi <- rbind(Fig3.DSB.Mi,aux.subset2)
  }
}


Fig3.ESB.Mi <- subset(ESB.Mi,select=c("SO2_BothOutcomes","NOX_BothOutcomes"))
aux.rownames <- paste0("BothOutcomes")
Fig3.ESB.Mi <- cbind(aux.rownames,Fig3.ESB.Mi)
colnames(Fig3.ESB.Mi) <- c("Data_Name","SO2","NOX")

#Showing values
cbind(Fig3.BEN.Mi[,1],round(rowSums(Fig3.BEN.Mi[,-1]),3))
cbind(Fig3.DSB.Mi[,1],round(rowSums(Fig3.DSB.Mi[,-1]),3))
cbind(Fig3.ESB.Mi[,1],round(rowSums(Fig3.ESB.Mi[,-1]),3))


###### 4.5.4. Results for Figure 4 ######
#a) DTOX -- ALL EFs: Benefits per mile, by outcome & urbanization
#b) DTOX -- ALL EFs: Health impacts of DSBs per mile, by pollutant & urbanization 
#c) ESBs: Health impacts of ESBs per mile, by pollutant
aux.ef <- c("NEI","MY2005","MY2010","MY2020")
aux.oc <- c("Deaths","Asthma")

for (i in 1:length(aux.ef))
{
  aux.subset <- paste("SUMDTOX",aux.ef[i],aux.oc,sep="_")
  aux.subset2 <- subset(DSB.Mi.URB, select=aux.subset)
  aux.subset2 <- aux.subset2[-7,]
  aux.d <- grep("Deaths",colnames(aux.subset2))
  aux.subset2[,aux.d] <- aux.subset2[,aux.d] - ESB.Mi$SUM_Deaths
  aux.a <- grep("Asthma",colnames(aux.subset2))
  aux.subset2[,aux.a] <- aux.subset2[,aux.a] - ESB.Mi$SUM_Asthma
  aux.names <- paste("SUMDTOX",aux.ef[i],Urb.Categories,sep="_") 
  aux.subset2 <- cbind(aux.names,aux.subset2)
  colnames(aux.subset2) <- c("Data_Name",aux.oc)
  if (i==1)
  {
    Fig4.BEN.Mi <- aux.subset2
  } else {
    Fig4.BEN.Mi <- rbind(Fig4.BEN.Mi,aux.subset2)
  }
}

aux.ef <- c("NEI","MY2005","MY2010","MY2020")
aux.pol <- c("PMDTOX","SO2","NOX","NH3","VOC")

for (i in 1:length(aux.ef))
{
  
  aux.subset <- paste(aux.pol,aux.ef[i],"BothOutcomes",sep="_")
  aux.subset2 <- subset(DSB.Mi.URB, select=aux.subset)
  aux.subset2 <- aux.subset2[-7,]
  aux.names <- paste(aux.ef[i],"BothOutcomes",Urb.Categories,sep="_") 
  aux.subset2 <- cbind(aux.names,aux.subset2)
  colnames(aux.subset2) <- c("Data_Name",aux.pol)
  if (i==1)
  {
    Fig4.DSB.Mi <- aux.subset2
  } else {
    Fig4.DSB.Mi <- rbind(Fig4.DSB.Mi,aux.subset2)
  }
}


Fig4.ESB.Mi <- subset(ESB.Mi,select=c("SO2_BothOutcomes","NOX_BothOutcomes"))
aux.rownames <- paste0("BothOutcomes")
Fig4.ESB.Mi <- cbind(aux.rownames,Fig4.ESB.Mi)
colnames(Fig4.ESB.Mi) <- c("Data_Name","SO2","NOX")

#Showing values
cbind(Fig4.BEN.Mi[,1],round(rowSums(Fig4.BEN.Mi[,-1]),3))
cbind(Fig4.DSB.Mi[,1],round(rowSums(Fig4.DSB.Mi[,-1]),3))
cbind(Fig4.ESB.Mi[,1],round(rowSums(Fig4.ESB.Mi[,-1]),3))


###### 4.5.5. Results for Figure S1 ######
#a) ALL EFs: Attrib deaths of DSBs per 1,000 buses, by pollutant & urbanization
#b) ESBs: Attrib deaths of ESBs per 1,000 buses, by pollutant
aux.ef <- c("NEI","MY2005","MY2010","MY2020")
aux.pol <- c("PM25","SO2","NOX","NH3","VOC")

for (i in 1:length(aux.ef))
{
  
  aux.subset <- paste(aux.pol,aux.ef[i],"AttribD",sep="_")
  aux.subset2 <- subset(DSB.BUS.URB, select=aux.subset)
  aux.subset2 <- aux.subset2[-7,]
  aux.subset2 <- 1000 * aux.subset2 #per 1,000 buses
  aux.names <- paste(aux.ef[i],"AttribD",Urb.Categories,sep="_") 
  aux.subset2 <- cbind(aux.names,aux.subset2)
  colnames(aux.subset2) <- c("Data_Name",aux.pol)
  if (i==1)
  {
    FigS1.DSB.BUS <- aux.subset2
  } else {
    FigS1.DSB.BUS <- rbind(FigS1.DSB.BUS,aux.subset2)
  }
}



FigS1.ESB.BUS <- subset(ESB.BUS,select=c("SO2_AttribD","NOX_AttribD"))
FigS1.ESB.BUS <- 1000 * FigS1.ESB.BUS #per 1,000 buses
aux.rownames <- paste0("AttribD")
FigS1.ESB.BUS <- cbind(aux.rownames,FigS1.ESB.BUS)
colnames(FigS1.ESB.BUS) <- c("Data_Name","SO2","NOX")

#Showing values
cbind(FigS1.DSB.BUS[,1],round(rowSums(FigS1.DSB.BUS[,-1]),3))
cbind(FigS1.ESB.BUS[,1],round(rowSums(FigS1.ESB.BUS[,-1]),3))

3.63*12.4*0.89

###### 4.5.6. Results for Figure S2 ######
#a) ALL EFs: Attrib new childhood asthma cases per 1,000 buses for DSBs, by pollutant & urbanization
#b) ESBs: Attrib new childhood asthma cases per 1,000 buses for ESBs, by pollutant
aux.ef <- c("NEI","MY2005","MY2010","MY2020")
aux.pol <- c("PM25","SO2","NOX","NH3","VOC")

for (i in 1:length(aux.ef))
{
  
  aux.subset <- paste(aux.pol,aux.ef[i],"AttribA",sep="_")
  aux.subset2 <- subset(DSB.BUS.URB, select=aux.subset)
  aux.subset2 <- aux.subset2[-7,]
  aux.subset2 <- 1000 * aux.subset2 #per 1,000 buses
  aux.names <- paste(aux.ef[i],"AttribA",Urb.Categories,sep="_") 
  aux.subset2 <- cbind(aux.names,aux.subset2)
  colnames(aux.subset2) <- c("Data_Name",aux.pol)
  if (i==1)
  {
    FigS2.DSB.BUS <- aux.subset2
  } else {
    FigS2.DSB.BUS <- rbind(FigS2.DSB.BUS,aux.subset2)
  }
}


FigS2.ESB.BUS <- subset(ESB.BUS,select=c("SO2_AttribA","NOX_AttribA"))
FigS2.ESB.BUS <- 1000 * FigS2.ESB.BUS #per 1,000 buses
aux.rownames <- paste0("AttribA")
FigS2.ESB.BUS <- cbind(aux.rownames,FigS2.ESB.BUS)
colnames(FigS2.ESB.BUS) <- c("Data_Name","SO2","NOX")

#Showing values
cbind(FigS2.DSB.BUS[,1],round(rowSums(FigS2.DSB.BUS[,-1]),3))
cbind(FigS2.ESB.BUS[,1],round(rowSums(FigS2.ESB.BUS[,-1]),3))


###### 4.5.7. Results for Figure S3 ######
#a) ALL EFs: Attrib deaths of DSBs per 100M miles, by pollutant & urbanization
#b) ESBs: Attrib deaths of ESBs per 100M miles, by pollutant
aux.ef <- c("NEI","MY2005","MY2010","MY2020")
aux.pol <- c("PM25","SO2","NOX","NH3","VOC")

for (i in 1:length(aux.ef))
{
  
  aux.subset <- paste(aux.pol,aux.ef[i],"AttribD",sep="_")
  aux.subset2 <- subset(DSB.Mi.URB, select=aux.subset)
  aux.subset2 <- aux.subset2[-7,]
  aux.subset2 <- 1e8 * aux.subset2 #per 100M miles
  aux.names <- paste(aux.ef[i],"AttribD",Urb.Categories,sep="_") 
  aux.subset2 <- cbind(aux.names,aux.subset2)
  colnames(aux.subset2) <- c("Data_Name",aux.pol)
  if (i==1)
  {
    FigS3.DSB.Mi <- aux.subset2
  } else {
    FigS3.DSB.Mi <- rbind(FigS3.DSB.Mi,aux.subset2)
  }
}


FigS3.ESB.Mi <- subset(ESB.Mi,select=c("SO2_AttribD","NOX_AttribD"))
FigS3.ESB.Mi <- 1e8 * FigS3.ESB.Mi #per 100M miles
aux.rownames <- paste0("AttribD")
FigS3.ESB.Mi <- cbind(aux.rownames,FigS3.ESB.Mi)
colnames(FigS3.ESB.Mi) <- c("Data_Name","SO2","NOX")

#Showing values
cbind(FigS3.DSB.Mi[,1],round(rowSums(FigS3.DSB.Mi[,-1]),3))
cbind(FigS3.ESB.Mi[,1],round(rowSums(FigS3.ESB.Mi[,-1]),3))

###### 4.5.8. Results for Figure S4 ######
#a) ALL EFs: Attrib new childhood asthma cases per 100M miles for DSBs, by pollutant & urbanization
#b) ESBs: Attrib new childhood asthma cases per 100M miles for ESBs, by pollutant
aux.ef <- c("NEI","MY2005","MY2010","MY2020")
aux.pol <- c("PM25","SO2","NOX","NH3","VOC")

for (i in 1:length(aux.ef))
{
  
  aux.subset <- paste(aux.pol,aux.ef[i],"AttribA",sep="_")
  aux.subset2 <- subset(DSB.Mi.URB, select=aux.subset)
  aux.subset2 <- aux.subset2[-7,]
  aux.subset2 <- 1e8 * aux.subset2 #per 100M miles
  aux.names <- paste(aux.ef[i],"AttribA",Urb.Categories,sep="_") 
  aux.subset2 <- cbind(aux.names,aux.subset2)
  colnames(aux.subset2) <- c("Data_Name",aux.pol)
  if (i==1)
  {
    FigS4.DSB.Mi <- aux.subset2
  } else {
    FigS4.DSB.Mi <- rbind(FigS4.DSB.Mi,aux.subset2)
  }
}


FigS4.ESB.Mi <- subset(ESB.Mi,select=c("SO2_AttribA","NOX_AttribA"))
FigS4.ESB.Mi <- 1e8 * FigS4.ESB.Mi #per 100M miles
aux.rownames <- paste0("AttribA")
FigS4.ESB.Mi <- cbind(aux.rownames,FigS4.ESB.Mi)
colnames(FigS4.ESB.Mi) <- c("Data_Name","SO2","NOX")

#Showing values
cbind(FigS4.DSB.Mi[,1],round(rowSums(FigS4.DSB.Mi[,-1]),3))
cbind(FigS4.ESB.Mi[,1],round(rowSums(FigS4.ESB.Mi[,-1]),3))

###### 4.5.9. Results for Figures S5-S8 ######
#Deaths & Asthma
aux.FigS5S8 <- cbind(Marg.Damages.GEMM_BasePM_2017_Mortality_2017,Marg.Damages.Tetreault_BasePM_2017_Asthma_2019[,c(16:20)])
#Sum
aux.sum <- aux.FigS5S8[,1:5] + aux.FigS5S8[,6:10]

FigS5S8.MD <- cbind(aux.FigS5S8,aux.sum)
aux.eg <- expand.grid(c("PM25","SO2","NOX","NH3","VOC"),c("Deaths","Asthma","BothOutcomes"))
colnames(FigS5S8.MD) <- paste("MD",aux.eg[,1],aux.eg[,2],sep="_")

#Showing values
summary(FigS5S8.MD)

#Fig 8 is weighted by Emissions
#Emissions are EFs * VMT 
#VMT is available at DSB.EFs.NEI$VMT (column 10)
#Emissions are available in columns 12-15 & 19 (for the Exhaust part of PM25)

#For example, for weighted means:
aux.emis.cols <- c(12,8:11)
aux.pollutant <- c("PM25","SO2","NOX","NH3","VOC")

for (aux.em in 1:5)
{
  aux.weight <- DSB.EFs.NEI[,aux.emis.cols[aux.em]] * DSB.EFs.NEI[,6] * 1e-6
  wt.means <- as.numeric(weighted_mean(FigS5S8.MD[,10+aux.em],w=aux.weight))
  wt.means <- round(wt.means,digits=0)
  print(paste(aux.pollutant[aux.em],"-- BothOutcomes --","Emissions-Weighted Mean = $",wt.means,"per metric ton",sep=" "))
  
  wt.means <- as.numeric(weighted_mean(FigS5S8.MD[,0+aux.em],w=aux.weight))
  wt.means <- round(wt.means,digits=0)
  print(paste(aux.pollutant[aux.em],"-- Deaths --","Emissions-Weighted Mean = $",wt.means,"per metric ton",sep=" "))
  
  wt.means <- as.numeric(weighted_mean(FigS5S8.MD[,5+aux.em],w=aux.weight))
  wt.means <- round(wt.means,digits=0)
  print(paste(aux.pollutant[aux.em],"-- Asthma --","Emissions-Weighted Mean = $",wt.means,"per metric ton",sep=" "))
}


###### 4.5.10. Results for Figure S9 ######
#a)NEI EFs: ESB benefits per mile, by outcome
#b)NEI EFs: DSB impacts per mile, by pollutant
FigS9.BEN.Mi <- subset(DSB.Mi, select=c(paste("SUM","NEI",c("Deaths","Asthma"),sep="_")))
FigS9.BEN.Mi$SUM_NEI_Deaths <- FigS9.BEN.Mi$SUM_NEI_Deaths - ESB.Mi$SUM_Deaths
FigS9.BEN.Mi$SUM_NEI_Asthma <- FigS9.BEN.Mi$SUM_NEI_Asthma - ESB.Mi$SUM_Asthma
names(FigS9.BEN.Mi)

FigS9.DSB.Mi <- subset(DSB.Mi, select=c(paste(c("PM25","SO2","NOX","NH3","VOC"),"NEI","BothOutcomes",sep="_")))
names(FigS9.DSB.Mi)

#Showing values
sum(rowSums(FigS9.BEN.Mi)*DSB.EFs.NEI$VMT)/sum(DSB.EFs.NEI$VMT) #VMT-weighted mean benefit
sum(rowSums(FigS9.DSB.Mi)*DSB.EFs.NEI$VMT)/sum(DSB.EFs.NEI$VMT) #VMT-weighted mean impact of DSBs


###### 4.5.11. Results for Figure S10 ######
#a) ALL EFs: Benefits per bus, by outcome & urbanization
#b) ALL EFs: Health impacts of DSBs per bus, by pollutant & urbanization 
#c) ESBs: Health impacts of ESBs per bus, by pollutant
aux.ef <- c("NEI","MY2005","MY2010","MY2020")
aux.oc <- c("Deaths","Asthma")

for (i in 1:length(aux.ef))
{
  aux.subset <- paste("SUM",aux.ef[i],aux.oc,sep="_")
  aux.subset2 <- subset(DSB.BUS.URB, select=aux.subset)
  aux.subset2 <- aux.subset2[-7,]
  aux.d <- grep("Deaths",colnames(aux.subset2))
  aux.subset2[,aux.d] <- aux.subset2[,aux.d] - ESB.BUS$SUM_Deaths
  aux.a <- grep("Asthma",colnames(aux.subset2))
  aux.subset2[,aux.a] <- aux.subset2[,aux.a] - ESB.BUS$SUM_Asthma
  aux.names <- paste("SUM",aux.ef[i],Urb.Categories,sep="_") 
  aux.subset2 <- cbind(aux.names,aux.subset2)
  colnames(aux.subset2) <- c("Data_Name",aux.oc)
  if (i==1)
  {
    FigS10.BEN.BUS <- aux.subset2
  } else {
    FigS10.BEN.BUS <- rbind(FigS10.BEN.BUS,aux.subset2)
  }
}

aux.ef <- c("NEI","MY2005","MY2010","MY2020")
aux.pol <- c("PM25","SO2","NOX","NH3","VOC")

for (i in 1:length(aux.ef))
{
  
  aux.subset <- paste(aux.pol,aux.ef[i],"BothOutcomes",sep="_")
  aux.subset2 <- subset(DSB.BUS.URB, select=aux.subset)
  aux.subset2 <- aux.subset2[-7,]
  aux.names <- paste(aux.ef[i],"BothOutcomes",Urb.Categories,sep="_") 
  aux.subset2 <- cbind(aux.names,aux.subset2)
  colnames(aux.subset2) <- c("Data_Name",aux.pol)
  if (i==1)
  {
    FigS10.DSB.BUS <- aux.subset2
  } else {
    FigS10.DSB.BUS <- rbind(FigS10.DSB.BUS,aux.subset2)
  }
}


FigS10.ESB.BUS <- subset(ESB.BUS,select=c("SO2_BothOutcomes","NOX_BothOutcomes"))
aux.rownames <- paste0("BothOutcomes")
FigS10.ESB.BUS <- cbind(aux.rownames,FigS10.ESB.BUS)
colnames(FigS10.ESB.BUS) <- c("Data_Name","SO2","NOX")

#Showing values
cbind(FigS10.BEN.BUS[,1],round(rowSums(FigS10.BEN.BUS[,-1]),-2))
cbind(FigS10.DSB.BUS[,1],round(rowSums(FigS10.DSB.BUS[,-1]),-2))
cbind(FigS10.ESB.BUS[,1],round(rowSums(FigS10.ESB.BUS[,-1]),-2))

###### 4.5.12. Results for Figure S11 ######
#a) DTOX -- ALL EFs: Benefits per bus, by outcome & urbanization
#b) DTOX -- ALL EFs: Health impacts of DSBs per bus, by pollutant & urbanization 
#c) ESBs: Health impacts of ESBs per bus, by pollutant
aux.ef <- c("NEI","MY2005","MY2010","MY2020")
aux.oc <- c("Deaths","Asthma")

for (i in 1:length(aux.ef))
{
  aux.subset <- paste("SUMDTOX",aux.ef[i],aux.oc,sep="_")
  aux.subset2 <- subset(DSB.BUS.URB, select=aux.subset)
  aux.subset2 <- aux.subset2[-7,]
  aux.d <- grep("Deaths",colnames(aux.subset2))
  aux.subset2[,aux.d] <- aux.subset2[,aux.d] - ESB.BUS$SUM_Deaths
  aux.a <- grep("Asthma",colnames(aux.subset2))
  aux.subset2[,aux.a] <- aux.subset2[,aux.a] - ESB.BUS$SUM_Asthma
  aux.names <- paste("SUMDTOX",aux.ef[i],Urb.Categories,sep="_") 
  aux.subset2 <- cbind(aux.names,aux.subset2)
  colnames(aux.subset2) <- c("Data_Name",aux.oc)
  if (i==1)
  {
    FigS11.BEN.BUS <- aux.subset2
  } else {
    FigS11.BEN.BUS <- rbind(FigS11.BEN.BUS,aux.subset2)
  }
}

aux.ef <- c("NEI","MY2005","MY2010","MY2020")
aux.pol <- c("PMDTOX","SO2","NOX","NH3","VOC")

for (i in 1:length(aux.ef))
{
  
  aux.subset <- paste(aux.pol,aux.ef[i],"BothOutcomes",sep="_")
  aux.subset2 <- subset(DSB.BUS.URB, select=aux.subset)
  aux.subset2 <- aux.subset2[-7,]
  aux.names <- paste(aux.ef[i],"BothOutcomes",Urb.Categories,sep="_") 
  aux.subset2 <- cbind(aux.names,aux.subset2)
  colnames(aux.subset2) <- c("Data_Name",aux.pol)
  if (i==1)
  {
    FigS11.DSB.BUS <- aux.subset2
  } else {
    FigS11.DSB.BUS <- rbind(FigS11.DSB.BUS,aux.subset2)
  }
}


FigS11.ESB.BUS <- subset(ESB.BUS,select=c("SO2_BothOutcomes","NOX_BothOutcomes"))
aux.rownames <- paste0("BothOutcomes")
FigS11.ESB.BUS <- cbind(aux.rownames,FigS11.ESB.BUS)
colnames(FigS11.ESB.BUS) <- c("Data_Name","SO2","NOX")

#Showing values
cbind(FigS11.BEN.BUS[,1],round(rowSums(FigS11.BEN.BUS[,-1]),-2))
cbind(FigS11.DSB.BUS[,1],round(rowSums(FigS11.DSB.BUS[,-1]),-2))
cbind(FigS11.ESB.BUS[,1],round(rowSums(FigS11.ESB.BUS[,-1]),-2))

### Check AttribD
cbind(FigS1.ESB.BUS[,1],round(rowSums(FigS1.ESB.BUS[,-1]),3))
cbind(FigS1.DSB.BUS[,1],round(rowSums(FigS1.DSB.BUS[,-1]),3))

sum(FigS1.DSB.BUS[7,-1]) - sum(FigS1.ESB.BUS[,-1])
sum(FigS2.DSB.BUS[7,-1]) - sum(FigS2.ESB.BUS[,-1])


###### 5. Plotting Figures 1, 3, 4, S1, S2, S3, S4, S10, S11 ######
#Creating labels
list.xlim <- list()
list.xaxs.at <- list()
list.xaxs.lab <- list()

#9 figures:
# Figures 1, 3, 4, S1, S2, S3, S4, S10, S11 
summary(rowSums(FigS1.DSB.BUS[,-1]))
summary(rowSums(FigS2.DSB.BUS[,-1]))

#x-axis: limits
list.xlim[[1]] <- 105000
list.xlim[[2]] <- 1.50
list.xlim[[3]] <- 5.25
list.xlim[[4]] <- 24
list.xlim[[5]] <- 44
list.xlim[[6]] <- 12.5
list.xlim[[7]] <- 23
list.xlim[[8]] <- 240000
list.xlim[[9]] <- 810000


#x-axis: placement of labels\
list.xaxs.at[[1]] <- 10000*0:10
list.xaxs.at[[2]] <- 0.25*0:6
list.xaxs.at[[3]] <- 0.5*0:10
list.xaxs.at[[4]] <- 2*0:11
list.xaxs.at[[5]] <- 5*0:8
list.xaxs.at[[6]] <- 1*0:12
list.xaxs.at[[7]] <- 2*0:11
list.xaxs.at[[8]] <- 20000*0:11
list.xaxs.at[[9]] <- 100000*0:8



#x-axis: labels
list.xaxs.lab[[1]] <- c("0",paste0("$",10*1:10,"k"))
list.xaxs.lab[[2]] <- paste0("$",format(0.25*0:6,nsmall=2))
list.xaxs.lab[[3]] <- paste0("$",format(0.5*0:10,nsmall=2))
list.xaxs.lab[[4]] <- 2*0:11
list.xaxs.lab[[5]] <- 5*0:8
list.xaxs.lab[[6]] <- 1*0:12
list.xaxs.lab[[7]] <- 2*0:11
list.xaxs.lab[[8]] <- c("0",paste0("$",20*1:11,"k"))
list.xaxs.lab[[9]] <- c("0",paste0("$",100*1:8,"k"))

#how many digits to round values to
list.round <- c(-2,3,3,2,2,2,2,-2,-2)


#titles -- 1st line
list.title.pt1 <- c(expression(bold("School Bus Health Impacts per bus and Benefits of Electrification: Fleet in 2017")),
                    expression(bold("School Bus Health Impacts per mile and Benefits of Electrification")),
                    expression(bold("School Bus Health Impacts per mile and Benefits of Electrification:")),
                    expression(bold(paste("School Bus: PM"["2.5"],"-attributable deaths per bus"))),
                    expression(bold(paste("School Bus: PM"["2.5"],"-attributable new childhood asthma cases per bus"))),
                    expression(bold(paste("School Bus: PM"["2.5"],"-attributable mortality per mile"))),
                    expression(bold(paste("School Bus: PM"["2.5"],"-attributable new childhood asthma cases per mile"))),
                    expression(bold("School Bus Health Impacts per bus and Benefits of Electrification")),
                    expression(bold("School Bus Health Impacts per bus:")))


#titles -- 2nd line
list.title.pt2 <- c(expression(bold("")),
                    expression(bold("")),
                    expression(bold("sensitivity analysis with carbonaceous particles 5 times more toxic by mass")),
                    expression(bold(paste(""))),
                    expression(bold(paste(""))),
                    expression(bold(paste(""))),
                    expression(bold(paste(""))),
                    expression(bold("")),
                    expression(bold("sensitivity analysis with carbonaceous particles 5 times more toxic by mass")))
                    
#x-axis: subtitle

list.xlab <- c(expression(bold("2022 USD per bus")),
               expression(bold("2022 USD per mile")),
               expression(bold("2022 USD per mile")),
               expression(bold(paste("PM"["2.5"],"-attributable deaths per 1,000 buses"))),
               expression(bold(paste("PM"["2.5"],"-attributable new childhood asthma cases per 1,000 buses"))),
               expression(bold(paste("PM"["2.5"],"-attributable deaths per 100 million miles driven"))),
               expression(bold(paste("PM"["2.5"],"-attributable new childhood asthma cases per 100 million miles driven"))),
               expression(bold("2022 USD per bus")),
               expression(bold("2022 USD per bus")))



#whether to display $ sign
list.k2 <- c("$",
             "$",
             "$",
             "",
             "",
             "",
             "",
             "$",
             "$")


#Names of urbanization categories
NCHS.Names <- c("Large Central Metro","Large Fringe Metro",
                "Medium Metro","Small Metro",
                "Micropolitan","Non-core ('Rural')","U.S. Average")

#Pollutant colors
col.pols <- c("blue","red","green3","darkorange","yellow")

#y axis positions
aux.mwd <- 0.4
aux.mywd <- 1
aux.ylim <- c(-3,(7*4)+3*aux.mywd)
y.gap <- 1
aux.ylim.f2 <- aux.ylim + c(0,y.gap+(aux.ylim[2]-aux.ylim[1])+y.gap)
f1.ygap1 <- 1.5
f1.ygap2 <- 1
aux.ylim.f1 <- c(0,1+1+f1.ygap1+7+f1.ygap1+1+f1.ygap1+7+f1.ygap2)



#File names
aux.fig.names <- c("1",
                   "3",
                   "4",
                   "S1",
                   "S2",
                   "S3",
                   "S4",
                   "S10",
                   "S11")

FigVersion <- "031124"
FigPath <- "Outputs/"

DSB.BUS.URB$SUM_NEI_BothOutcomes
rowSums(Fig1.BEN.BUS[,-1])[7] + BEN.BUS.GHG


list.fg.n <- paste0("Fig_",aux.fig.names,"_",FigVersion,".pdf")
list.fg.n <- paste0(FigPath,list.fg.n)


Fig.BEN.Datasets <- c("Fig1.BEN.BUS","Fig3.BEN.Mi","Fig4.BEN.Mi",
                      "","","","",
                      "FigS10.BEN.BUS","FigS11.BEN.BUS")
Fig.DSB.Datasets <- c("Fig1.DSB.BUS","Fig3.DSB.Mi","Fig4.DSB.Mi",
                      "FigS1.DSB.BUS","FigS2.DSB.BUS","FigS3.DSB.Mi","FigS4.DSB.Mi",
                      "FigS10.DSB.BUS","FigS11.DSB.BUS")
Fig.ESB.Datasets <- c("Fig1.ESB.BUS","Fig3.ESB.Mi","Fig4.ESB.Mi",
                      "FigS1.ESB.BUS","FigS2.ESB.BUS","FigS3.ESB.Mi","FigS4.ESB.Mi",
                      "FigS10.ESB.BUS","FigS11.ESB.BUS")

# Fig1.ESB.BUS
# Fig3.ESB.Mi
# Fig4.ESB.Mi
# FigS1.ESB.BUS
# FigS2.ESB.BUS
# FigS3.ESB.Mi
# FigS4.ESB.Mi
# FigS10.ESB.BUS
# FigS11.ESB.BUS

aux.old.i <- c(10,9,11,4,5,7,8,1,12)

# fg <- 9



for (fg in c(2:9))
{
  if(!(fg %in% c(2,3,8,9)))
  {
    pdf(file=list.fg.n[fg],width=12,height=8.4)
    aux.ylim.plot <- aux.ylim
    par(xaxs="i",yaxs="i",
        mar=c(3.5,11.5,3.25,1))
  } else {
    pdf(file=list.fg.n[fg],width=12,height=8.4*(diff(aux.ylim.f2)/diff(aux.ylim)))
    aux.ylim.plot <- aux.ylim.f2
    par(xaxs="i",yaxs="i",
        mar=c(3.5,13,3.25,1))
  }
  
  plot(x=-100,y=-100,
       xlab="",ylab="",
       xaxt="n",yaxt="n",
       xlim=c(0,list.xlim[[fg]]),ylim=aux.ylim.plot)
  
  aux.y1 <- aux.ylim[2] - (rep(c(1:7),4) + rep(((7+aux.mywd)*0:3),7))
  aux.y1 <- c((aux.ylim[1]+1),aux.y1)
  
  aux.y2 <- aux.ylim[2] - (4 + ((7+aux.mywd)*0:3))
  mtext(text=c("DSB: Fleet in 2017", "DSB: MY 2005","DSB: MY 2010","DSB: MY 2020"),at=aux.y2,
        side=2,line=10,font=2,cex=1.2)
  
  mtext(text=list.xlab[fg],at=mean(c(0,list.xlim[[fg]])),
        side=1,line=2.25,font=2,cex.axis=1.2)
  
  axis(side=2,at=c(aux.y1[-1]),
       labels=c(rep(NCHS.Names,4)),las=2,cex.axis=1.2)
  
  axis(side=2,at=c(aux.y1[1]),
       labels=c("Electric School Bus"),las=2, font=2, cex.axis=1.2)
  
  axis(side=1,at=list.xaxs.at[[fg]],
       labels=list.xaxs.lab[[fg]])
  
  if(fg %in% c(3,9))
  {
    mtext(text=list.title.pt2[fg],
          side=3,line=0.3,font=2,cex=1.2, at=mean(c(0,list.xlim[[fg]])))
    
    mtext(text=list.title.pt1[fg],
          side=3,line=1.6,font=2,cex=1.2, at=mean(c(0,list.xlim[[fg]])))
  } else {
    mtext(text=list.title.pt1[fg],
          side=3,line=0.7,font=2,cex=1.4, at=mean(c(0,list.xlim[[fg]])))
  }
  
  if(fg %in% c(1:3,8,9))
  {
    BEN.Data <- get(Fig.BEN.Datasets[fg])
    BEN.Data[,1] <- 0


  }
  
  
  #Plotting DSB Data
  DSB.Data <- get(Fig.DSB.Datasets[fg])
  DSB.Data[,1] <- 0
  DSB.Data[,2:6] <- DSB.Data[,c(2,4,6,5,3)]
  
  aux.rec.y <- c(30:24,22:16,14:8,6:0)
  for(j in 1:5)
  {
    rect(xleft=rowSums(as.matrix(DSB.Data[,1:j])),
         xright=rowSums(as.matrix(DSB.Data[,1:(j+1)])),
         ybottom=aux.rec.y-aux.mwd,
         ytop=aux.rec.y+aux.mwd,
         col=col.pols[j])
  }

  aux.text <- format(as.numeric(round(rowSums(DSB.Data[1,]),digits=list.round[fg])),nsmall=max(0,list.round[fg]),big.mark = ",")
  for(tx in 2:dim(DSB.Data)[1])
  {
    aux.text1 <- format(as.numeric(round(rowSums(DSB.Data[tx,]),digits=list.round[fg])),nsmall=max(0,list.round[fg]),big.mark = ",")
    aux.text <- c(aux.text,aux.text1)
  }
  
  if(fg %in% c(2,3))
  {
    text.d0 <- aux.text
    for (aux.d in 1:length(aux.text))
    {
      text.d1 <- as.numeric(text.d0[aux.d])
      if(text.d1<0.095)
      {
        text.d1 <- format(text.d1,digits=2,nsmall=3)
      } else {
        text.d1 <- format(text.d1,digits=2,nsmall=2)
      }
      
      aux.text[aux.d] <- text.d1 
    }
  }
  aux.text.pos <- as.numeric(rowSums(DSB.Data)) + (0.01*list.xlim[[fg]])
  for(tx in 1:length(aux.text))
  {
    text(paste0(list.k2[fg],aux.text[tx]),y=aux.rec.y[tx],x=aux.text.pos[tx],adj=c(0,0.5),cex=1.2)
  }

  
  #Plotting ESB Data
  ESB.Data <- get(Fig.ESB.Datasets[fg])
  ESB.Data[,1] <- 0
  ESB.Data[,2:3] <- ESB.Data[,3:2]

  col.pols2 <- col.pols[c(2,5)]

  for(k in 1:2)
  {
    rect(xleft=rowSums(as.matrix(ESB.Data[c(1:k)])),
         xright=rowSums(as.matrix(ESB.Data[c(1:(k+1))])),
         ybottom=aux.ylim.plot[1]+1-aux.mwd,
         ytop=aux.ylim.plot[1]+1+aux.mwd,
         col=col.pols2[k])
    
  }

  aux.text <- format(as.numeric(round(rowSums(as.matrix(ESB.Data)),digits=list.round[fg])),nsmall=max(0,list.round[fg]),big.mark = ",")

  if(fg %in% c(2,3))
  {
    text.d0 <- aux.text
    for (aux.d in 1:length(aux.text))
    {
      text.d1 <- as.numeric(text.d0[aux.d])
      if(text.d1<0.095)
      {
        text.d1 <- format(text.d1,digits=2,nsmall=3)
      } else {
        text.d1 <- format(text.d1,digits=2,nsmall=2)
      }
      
      aux.text[aux.d] <- text.d1 
    }
  }
  aux.text.pos <- as.numeric(rowSums(as.matrix(ESB.Data))) + (0.01*list.xlim[[fg]])
  text(paste0(list.k2[fg],aux.text),y=aux.ylim.plot[1]+1,x=aux.text.pos,adj=c(0,0.5),cex=1.2)

  if (!(fg %in% c(2,3,8,9)))
  {
    legend(x=list.xlim[[fg]]*0.99,
           y=(aux.ylim[2]-1),
           xjust=1,
           yjust=1,
           legend=c(expression(paste("PM"["2.5"])),
                    expression(paste("NO"["x"])),
                    expression(paste("VOC")),
                    expression(paste("NH"["3"])),
                    expression(paste("SO"["2"]))),
           pch=22,pt.bg=col.pols,col="black",bty="n",cex=1.3,pt.cex=1.3)
  }
  
  if(fg %in% c(2,3,8,9))
  {
    
    legend(x=list.xlim[[fg]]*0.99,
           y=(aux.ylim[2]+y.gap-1),
           xjust=1,
           yjust=1,
           legend=c(expression(paste("PM"["2.5"])),
                    expression(paste("NO"["x"])),
                    expression(paste("VOC")),
                    expression(paste("NH"["3"])),
                    expression(paste("SO"["2"]))),
           pch=22,pt.bg=col.pols,col="black",bty="n",cex=1.3,pt.cex=1.3)
    
    if (fg %in% c(2,3))
    {
      text("(b) Health Impacts Per Mile Driven",
           x=0.01*list.xlim[[fg]],
           y=aux.ylim[2],
           adj=c(0,0),
           cex=1.3, font=2)
      
      text("(a) ESB Benefits Per Mile Driven",
           x=0.01*list.xlim[[fg]],
           y=aux.ylim.f2[2]-y.gap,
           adj=c(0,0),
           cex=1.3, font=2)
    } else {
      text("(b) Health Impacts Per Bus",
           x=0.01*list.xlim[[fg]],
           y=aux.ylim[2],
           adj=c(0,0),
           cex=1.3, font=2)
      
      text("(a) ESB Benefits Per Bus",
           x=0.01*list.xlim[[fg]],
           y=aux.ylim.f2[2]-y.gap,
           adj=c(0,0),
           cex=1.3, font=2)
    }
    
    aux.y1t <- aux.ylim.f2[2] - y.gap - (rep(c(1:7),4) + rep(((7+aux.mywd)*0:3),7))
    aux.y1t <- c((aux.ylim[2]+y.gap+1),aux.y1t)
    
    aux.y2t <- aux.ylim.f2[2] - y.gap - (4 + ((7+aux.mywd)*0:3))
    mtext(text=c("DSB: Fleet in 2017", "DSB: MY 2005","DSB: MY 2010","DSB: MY 2020"),at=aux.y2t,
          side=2,line=10,font=2,cex=1.2)
    
    mtext(text=c("Health Benefits of Replacing Different DSBs with ESBs"),at=mean(range(aux.y1t)),
          side=2,line=11.5,font=2,cex=1.3)
    
    axis(side=2,at=c(aux.y1t[-1]),
         labels=c(rep(NCHS.Names,4)),las=2,cex.axis=1.2)
    
    axis(side=2,at=c(aux.y1t[1]),
         labels=c("Climate Benefits"),las=2, font=2, cex.axis=1.3)
    
    aux.y1tb <- aux.y1t[-1]
    aux.y1tb <- aux.y1tb[order(aux.y1tb,decreasing=TRUE)]
    
    col.outcomes <- c("turquoise","magenta","slateblue")

    for (j in 1:2)
    {
      rect(xleft=rowSums(as.matrix(BEN.Data[,1:j])),
           xright=rowSums(as.matrix(BEN.Data[,1:(j+1)])),
           ybottom=aux.y1tb-aux.mwd,
           ytop=aux.y1tb+aux.mwd,
           col=col.outcomes[j])
    }
    
    aux.text <- format(as.numeric(round(rowSums(BEN.Data[1,]),digits=list.round[fg])),nsmall=max(0,list.round[fg]),big.mark = ",")
    for(tx in 2:dim(BEN.Data)[1])
    {
      aux.text1 <- format(as.numeric(round(rowSums(BEN.Data[tx,]),digits=list.round[fg])),nsmall=max(0,list.round[fg]),big.mark = ",")
      aux.text <- c(aux.text,aux.text1)
    }
    
    if(fg %in% c(2,3))
    {
      text.d0 <- aux.text
      for (aux.d in 1:length(aux.text))
      {
        text.d1 <- as.numeric(text.d0[aux.d])
        if(text.d1<0.095)
        {
          text.d1 <- format(text.d1,digits=2,nsmall=3)
        } else {
          text.d1 <- format(text.d1,digits=2,nsmall=2)
        }
        
        aux.text[aux.d] <- text.d1 
      }
    }
    
    aux.text.pos <- as.numeric(rowSums(BEN.Data)) + (0.01*list.xlim[[fg]])

    for(tx in 1:length(aux.text))
    {
      text(paste0(list.k2[fg],aux.text[tx]),y=aux.y1tb[tx],x=aux.text.pos[tx],adj=c(0,0.5),cex=1.2)
    }
    
    aux.GHG2 <- BEN.BUS.GHG
    
    if(fg %in% c(2,3))
    {
      aux.GHG2 <- BEN.Mi.GHG
    }
    
    rect(xleft=0,
         xright=aux.GHG2,
         ybottom=aux.y1t[1]-aux.mwd,
         ytop=aux.y1t[1]+aux.mwd,
         col=col.outcomes[3])
    
    aux.text <- format(round(as.numeric(aux.GHG2),digits=list.round[fg]),nsmall=max(0,list.round[fg]),big.mark = ",")
    
    if(fg %in% c(2,3))
    {
      text.d0 <- aux.text
      for (aux.d in 1:length(aux.text))
      {
        text.d1 <- as.numeric(text.d0[aux.d])
        if(text.d1<0.095)
        {
          text.d1 <- format(text.d1,digits=2,nsmall=3)
        } else {
          text.d1 <- format(text.d1,digits=2,nsmall=2)
        }
        
        aux.text[aux.d] <- text.d1 
      }
    }
    
    aux.text.pos <- as.numeric(aux.GHG2) + (0.01*list.xlim[[fg]])
    text(paste0(list.k2[fg],aux.text),y=aux.y1t[1],x=aux.text.pos,adj=c(0,0.5),cex=1.2)
    
    legend(x=list.xlim[[fg]]*0.99,
           y=(aux.ylim.f2[2]-1),
           xjust=1,
           yjust=1,
           legend=c("Mortality",
                    "Childhood Asthma",
                    "Climate Change"),
           pch=22,pt.bg=col.outcomes,col="black",bty="n",cex=1.3,pt.cex=1.3)
    
  }
  
  dev.off()
}

f1.ygap1 <- 1.5
f1.ygap2 <- 2

aux.ylim.f1 <- c(0,1+1+f1.ygap1+6+f1.ygap2+1+f1.ygap1+6+f1.ygap2*0.35)

range(aux.ylim.f1)
diff(aux.ylim)

for (fg in c(1))
{
  
  pdf(file=list.fg.n[fg],width=0.85*12,height=8.4*1.75*0.85*(diff(aux.ylim.f1)/diff(aux.ylim)))
  aux.ylim.plot <- aux.ylim.f1
  par(xaxs="i",yaxs="i",
      mar=c(3.5,11.5,3.25,1))
  
  plot(x=-100,y=-100,
       xlab="",ylab="",
       xaxt="n",yaxt="n",
       xlim=c(0,list.xlim[[fg]]),ylim=aux.ylim.plot)
  
  
  aux.y1 <- aux.ylim.f1[1] + c(1,c(1+f1.ygap1+0:6),c(1+f1.ygap1+6+f1.ygap2+1),c(1+f1.ygap1+6+f1.ygap2+1+f1.ygap1+0:6))

  mtext(text=c(rep("DSB: Fleet in 2017",2)),at=c(mean(aux.y1[c(10,16)]),mean(aux.y1[c(2,8)])),
        side=2,line=10,font=2,cex=1.2)

  mtext(text=list.xlab[fg],at=mean(c(0,list.xlim[[fg]])),
        side=1,line=2.25,font=2,cex=1.2)
  
  axis(side=2,at=c(aux.y1[-c(1,9)]),
       labels=c(rep(rev(NCHS.Names),2)),las=2,cex.axis=1.2)
  
  axis(side=2,at=c(aux.y1[1]),
       labels=c("Electric School Bus"),las=2, font=2, cex.axis=1.2)
  
  axis(side=2,at=c(aux.y1[9]),
       labels=c("Climate Benefits"),las=2, font=2, cex.axis=1.2)
  
  axis(side=1,at=list.xaxs.at[[fg]],
       labels=list.xaxs.lab[[fg]], cex.axis=1.2)
  
  mtext(text=list.title.pt2[fg],
        side=3,line=0.3,font=2,cex=1.2, at=mean(c(0,list.xlim[[fg]])))
  
  mtext(text=list.title.pt1[fg],
        side=3,line=1.6,font=2,cex=1.2, at=mean(c(0,list.xlim[[fg]])))
  
  #DSB
  DSB.Data <- get(Fig.DSB.Datasets[fg])
  DSB.Data[,1] <- 0
  DSB.Data[,2:6] <- DSB.Data[,c(2,4,6,5,3)]
  
  for (j in 1:5)
  {
    rect(xleft=rowSums(as.matrix(DSB.Data[,1:j])),
         xright=rowSums(as.matrix(DSB.Data[,1:(j+1)])),
         ybottom=aux.y1[8:2]-aux.mwd,
         ytop=aux.y1[8:2]+aux.mwd,
         col=col.pols[j])
  }
  
  aux.text <- format(as.numeric(round(rowSums(DSB.Data[1,]),digits=list.round[fg])),nsmall=max(0,list.round[fg]),big.mark = ",")
  for(tx in 2:dim(DSB.Data)[1])
  {
    aux.text1 <- format(as.numeric(round(rowSums(DSB.Data[tx,]),digits=list.round[fg])),nsmall=max(0,list.round[fg]),big.mark = ",")
    aux.text <- c(aux.text,aux.text1)
  }
  
  aux.text.pos <- as.numeric(rowSums(DSB.Data)) + (0.01*list.xlim[[fg]])
  for(tx in 1:length(aux.text))
  {
    text(paste0(list.k2[fg],aux.text[tx]),y=aux.y1[9-tx],x=aux.text.pos[tx],adj=c(0,0.5),cex=1.2)
  }
  
  #ESB
  ESB.Data <- get(Fig.ESB.Datasets[fg])
  ESB.Data[,1] <- 0
  ESB.Data[,2:3] <- ESB.Data[,3:2]
  
  col.pols2 <- col.pols[c(2,5)]
  
  for(k in 1:2)
  {
    rect(xleft=rowSums(as.matrix(ESB.Data[c(1:k)])),
         xright=rowSums(as.matrix(ESB.Data[c(1:(k+1))])),
         ybottom=aux.ylim.plot[1]+1-aux.mwd,
         ytop=aux.ylim.plot[1]+1+aux.mwd,
         col=col.pols2[k])
    
  }
  
  aux.text <- format(as.numeric(round(rowSums(as.matrix(ESB.Data)),digits=list.round[fg])),nsmall=max(0,list.round[fg]),big.mark = ",")
  
  aux.text.pos <- as.numeric(rowSums(as.matrix(ESB.Data))) + (0.01*list.xlim[[fg]])
  text(paste0(list.k2[fg],aux.text),y=aux.ylim.plot[1]+1,x=aux.text.pos,adj=c(0,0.5),cex=1.2)
  
  col.outcomes <- c("turquoise","magenta","slateblue")
  
  BEN.Data <- get(Fig.BEN.Datasets[fg])
  BEN.Data[,1] <- 0
  
  for (j in 1:2)
  {
    rect(xleft=rowSums(as.matrix(BEN.Data[,1:j])),
         xright=rowSums(as.matrix(BEN.Data[,1:(j+1)])),
         ybottom=aux.y1[16:10]-aux.mwd,
         ytop=aux.y1[16:10]+aux.mwd,
         col=col.outcomes[j])
  }
  
  aux.text <- format(as.numeric(round(rowSums(BEN.Data[1,]),digits=list.round[fg])),nsmall=max(0,list.round[fg]),big.mark = ",")
  for(tx in 2:dim(BEN.Data)[1])
  {
    aux.text1 <- format(as.numeric(round(rowSums(BEN.Data[tx,]),digits=list.round[fg])),nsmall=max(0,list.round[fg]),big.mark = ",")
    aux.text <- c(aux.text,aux.text1)
  }
  
  aux.text.pos <- as.numeric(rowSums(BEN.Data)) + (0.01*list.xlim[[fg]])

  for(tx in 1:length(aux.text))
  {
    text(paste0(list.k2[fg],aux.text[tx]),y=aux.y1[17-tx],x=aux.text.pos[tx],adj=c(0,0.5),cex=1.2)
  }
  
  aux.GHG2 <- BEN.BUS.GHG
  
  rect(xleft=0,
       xright=aux.GHG2,
       ybottom=aux.y1[9]-aux.mwd,
       ytop=aux.y1[9]+aux.mwd,
       col=col.outcomes[3])
  
  aux.text <- format(round(as.numeric(aux.GHG2),digits=list.round[fg]),nsmall=max(0,list.round[fg]),big.mark = ",")
  aux.text.pos <- as.numeric(aux.GHG2) + (0.01*list.xlim[[fg]])
  text(paste0(list.k2[fg],aux.text),y=aux.y1[9],x=aux.text.pos,adj=c(0,0.5),cex=1.2)
  
  legend(x=list.xlim[[fg]]*0.99,
         y=(aux.y1[1]-1),
         xjust=1,
         yjust=0,
         legend=c(expression(paste("PM"["2.5"])),
                  expression(paste("NO"["x"])),
                  expression(paste("VOC")),
                  expression(paste("NH"["3"])),
                  expression(paste("SO"["2"]))),
         pch=22,pt.bg=col.pols,col="black",bty="n",cex=1.3,pt.cex=1.3)
  
  
  legend(x=list.xlim[[fg]]*0.99,
         y=(aux.y1[9]-1),
         xjust=1,
         yjust=0,
         legend=c("Mortality",
                  "Childhood Asthma",
                  "Climate Change"),
         pch=22,pt.bg=col.outcomes,col="black",bty="n",cex=1.3,pt.cex=1.3)
  
  text("(b) Health Impacts Per Bus",
       x=0.01*list.xlim[[fg]],
       y=aux.y1[8]+0.8,
       adj=c(0,0),
       cex=1.3, font=2)
  
  text("(a) ESB Benefits Per Bus",
       x=0.01*list.xlim[[fg]],
       y=aux.y1[16]+0.8,
       adj=c(0,0),
       cex=1.3, font=2)
  
  dev.off()
}

###### 6. Plotting Figures S5, S6, S7, S8 ######

Plot.titles <- paste0("CDFs of marginal damages for each species -- ",
                      c("Mortality + Childhood Asthma","Mortality","Childhood Asthma"))
Plot.titles2 <- c("[weighted by total 2017 School Bus Emissions of each species]")

Plot.titles3 <- paste0("Distribution of marginal damages for each species by county -- ",
                       c("Mortality + Childhood Asthma","Mortality","Childhood Asthma"))
Plot.titles4 <- c("[unweighted, for the 3,108 counties in the contiguous U.S.]")

aux.fig.names2 <- c("S5","S6","S7","S8")

Plot.names <- paste0("Fig_",aux.fig.names2,"_",FigVersion,".pdf")
Plot.names <- paste0(FigPath,Plot.names)

aux.col.pol <- col.pols[c(1,5,2,4,3)]
aux.col.outcomes <- c("royalblue","turquoise","magenta")

#Figures S5-S7

for (crf in c(1,2,3))
{
  if(crf==1)
  {
    aux.marg.crf <- FigS5S8.MD[,11:15]
  }
  
  if(crf==2)
  {
    aux.marg.crf <- FigS5S8.MD[,1:5]
  }
  
  if(crf==3)
  {
    aux.marg.crf <- FigS5S8.MD[,6:10]
  }

  aux.main.text <- c(expression(bold(paste("Primary PM"["2.5"],"",sep=""))),
                     expression(bold(paste("SO"["2"],"",sep=""))),
                     expression(bold(paste("NO"["x"],"",sep=""))),
                     expression(bold(paste("NH"["3"],"",sep=""))),
                     expression(bold(paste("VOC", "",sep=""))))
  
  pdf(Plot.names[crf],height=10,width=6.5)
  par(mfcol=c(3,2),mgp=c(2,0.5,0),mar=c(3,3,2,1),
      cex.axis=0.8,cex.lab=0.8,oma=c(0.1,0.1,3.5,0.1),
      xaxs="i", yaxs="i")
  
  for (aux.pol in c(1,5,2,3,4))
  {
    aux.plot <- as.data.frame(matrix(rep(0,3108*2),ncol=2))
    names(aux.plot) <- c("MD","EMIS")
    aux.plot$MD <- log10(aux.marg.crf[,aux.pol])
    
    10^c(0.25*4:29)
    
    hist(aux.plot$MD,breaks=c(0.25*4:28),
         xlim=log10(c(10,1e7)),ylim=c(0,1700),
         xaxt="n",yaxt="n",
         xlab=expression(bold("Marginal Damage [2022 USD per metric ton], log scale")),
         ylab=expression(bold("Number of Counties")),
         main=aux.main.text[aux.pol],col=aux.col.pol[aux.pol],border="black")
    mtext(side=3, outer=TRUE, Plot.titles3[crf], line=1.9, font=2, cex=0.85)
    mtext(side=3, outer=TRUE, Plot.titles4, line=0.6, font=2,cex=0.75)
    
    axis(side=1, at=c(1:7),labels=c("$10","$100","$1k","$10k","$100k","$1M","$10M"),font=2)
    aux.ylab <- 200*0:8
    axis(side=2, at=aux.ylab,
         labels=aux.ylab,
         las=2,font=2)
    
    aux.ylab <- c(0.025,0.1,0.25,0.5,0.75,0.9,0.975)
    aux.x1 <- 10^(c(min(aux.plot$MD),quantile(aux.plot$MD, probs=aux.ylab),max(aux.plot$MD)))
    
    aux.MD.legend <- rep("a",length(aux.x1))
    for (aux.leg in 1:length(aux.x1))
    {
      
      if(aux.x1[aux.leg]<9950000)
      {
        aux.MD.legend[aux.leg] <- paste0(as.character(format(round(aux.x1[aux.leg]/1e6,digits=2)),nsmall=2),"M")
      }
      
      if(aux.x1[aux.leg]<995000)
      {
        aux.MD.legend[aux.leg] <- paste0(as.character(format(round(aux.x1[aux.leg]/1e3,digits=0)),nsmall=0),"k")
      }
      
      if(aux.x1[aux.leg]<99500)
      {
        aux.MD.legend[aux.leg] <- paste0(as.character(format(round(aux.x1[aux.leg]/1e3,digits=1)),nsmall=1),"k")
      }
      
      if(aux.x1[aux.leg]<9950)
      {
        aux.MD.legend[aux.leg] <- paste0(as.character(format(round(aux.x1[aux.leg]/1e3,digits=2)),nsmall=2),"k")
      }
      
      if(aux.x1[aux.leg]<995)
      {
        aux.MD.legend[aux.leg] <- paste0(as.character(format(round(aux.x1[aux.leg],digits=0)),nsmall=0),"")
      }
      
      if(aux.x1[aux.leg]<99.5)
      {
        aux.MD.legend[aux.leg] <- paste0(as.character(format(round(aux.x1[aux.leg],digits=1)),nsmall=0),"")
      }
      
    }
    
    aux.MD.legend2 <- c(eval(bquote(expression(bold(paste("Percentiles:  ",sep=""))))),
                        eval(bquote(expression(bold(paste("Min.:  ",.(aux.MD.legend[1]),sep=""))))),
                        eval(bquote(expression(bold(paste("2.5"^"th",":  ",.(aux.MD.legend[2]),sep=""))))),
                        eval(bquote(expression(bold(paste("10"^"th",":  ",.(aux.MD.legend[3]),sep=""))))),
                        eval(bquote(expression(bold(paste("25"^"th",":  ",.(aux.MD.legend[4]),sep=""))))),
                        eval(bquote(expression(bold(paste("50"^"th",":  ",.(aux.MD.legend[5]),sep=""))))),
                        eval(bquote(expression(bold(paste("75"^"th",":  ",.(aux.MD.legend[6]),sep=""))))),
                        eval(bquote(expression(bold(paste("90"^"th",":  ",.(aux.MD.legend[7]),sep=""))))),
                        eval(bquote(expression(bold(paste("97.5"^"th",":  ",.(aux.MD.legend[8]),sep=""))))),
                        eval(bquote(expression(bold(paste("Max.:  ",.(aux.MD.legend[9]),sep=""))))))
    
    legend("topright",
           legend=aux.MD.legend2, bty="n", cex=0.7, 
           pch=NA,col=NA,text.font=2)
    box(which="plot",lty="solid")
  }
  dev.off()
}



#Figure S8

aux.main.text <- c(expression(bold(paste("Primary PM"["2.5"],"",sep=""))),
                   expression(bold(paste("SO"["2"],"",sep=""))),
                   expression(bold(paste("NO"["x"],"",sep=""))),
                   expression(bold(paste("NH"["3"],"",sep=""))),
                   expression(bold(paste("VOC", "",sep=""))))

emis.col <- c(12,8:11)
emis.col.vmt <- 6

plot(0)
pdf(Plot.names[4],height=10,width=6.5)
par(mfcol=c(3,2),mgp=c(2,0.5,0),mar=c(3,3,2,0.5),
    cex.axis=0.8,cex.lab=0.8,oma=c(0.1,0.1,3.5,0.1))

for (aux.pol in c(1,5,2,3,4))
{
  plot(NA,xlim=log10(c(10,1e7)),ylim=c(0,1),
       xaxt="n",yaxt="n",
       xlab=expression(bold("Marginal Damage [2022 USD per metric ton], log scale")),
       ylab=expression(bold("CDF")),
       main=aux.main.text[aux.pol])
  
  mtext(side=3, outer=TRUE, "CDFs of marginal damages for each species", line=1.9, font=2)
  mtext(side=3, outer=TRUE, Plot.titles2, line=0.6, font=2,cex=0.8)
  
  axis(side=1, at=c(1:7),labels=c("$10","$100","$1k","$10k","$100k","$1M","$10M"),font=2)
  aux.ylab <- c(0,0.1,0.25,0.5,0.75,0.9,1)
  axis(side=2, at=aux.ylab,
       labels=aux.ylab,
       las=2,font=2)
  
  for (crf in c(1,2,3))
  {
    
    if(crf==1)
    {
      aux.marg.crf <- FigS5S8.MD[,11:15]
    }
    
    if(crf==2)
    {
      aux.marg.crf <- FigS5S8.MD[,1:5]
    }
    
    if(crf==3)
    {
      aux.marg.crf <- FigS5S8.MD[,6:10]
    }
    
    aux.plot <- as.data.frame(matrix(rep(0,3108*2),ncol=2))
    names(aux.plot) <- c("MD","EMIS")
    aux.plot$MD <- log10(aux.marg.crf[,aux.pol])
    
    aux.plot$EMIS <- (DSB.EFs.NEI[,emis.col[aux.pol]] * DSB.EFs.NEI[,emis.col.vmt]) * (1/1e6)
    aux.plot <- aux.plot[order(aux.plot$MD, decreasing=FALSE),]
    
    aux.plot$EMIS.Left <- 0
    aux.plot$EMIS.Right <- 0
    aux.plot$EMIS.Right[1] <- aux.plot$EMIS[1]/sum(aux.plot$EMIS)
    
    for (i in 2:length(aux.plot$MD))
    {
      aux.plot$EMIS.Left[i] <- sum(aux.plot$EMIS[1:(i-1)])/sum(aux.plot$EMIS)
      aux.plot$EMIS.Right[i] <- sum(aux.plot$EMIS[1:(i)])/sum(aux.plot$EMIS)
    }
    
    points(x=aux.plot$MD,y=aux.plot$EMIS.Right, pch=20, col=aux.col.outcomes[crf], cex=0.3)
    segments(x0=aux.plot$MD[-length(aux.plot$MD)],
             x1=aux.plot$MD[-1],
             y0=aux.plot$EMIS.Right[-length(aux.plot$MD)],
             y1=aux.plot$EMIS.Right[-length(aux.plot$MD)],
             lwd=0.5,col=aux.col.outcomes[crf])
    
    
    aux.OnRoad.Mean <- sum(10^aux.plot$MD * aux.plot$EMIS)/sum(aux.plot$EMIS)
    aux.over <- which.max(10^aux.plot$MD>aux.OnRoad.Mean)
    aux.over2 <- (10^aux.plot$MD[aux.over] - aux.OnRoad.Mean)/(10^aux.plot$MD[aux.over] - 10^aux.plot$MD[aux.over-1])
    aux.y.OnRoad.Mean <- aux.plot$EMIS.Right[aux.over-1] + ((1-aux.over2) * (aux.plot$EMIS.Right[aux.over] - aux.plot$EMIS.Right[aux.over-1]))
    points(x=log10(aux.OnRoad.Mean),y=aux.y.OnRoad.Mean,col="black",pch=8,cex=1)
    
    print(aux.OnRoad.Mean)
  }
  
}

plot(NA,xlim=c(0,1),ylim=c(0,1),
     xaxt="n",yaxt="n",
     xlab="",
     ylab="",
     bty="n")

legend(c("Mortality + Asthma","Mortality","Asthma","CDF Means"),
       pch=c(15,15,15,8), col=c(aux.col.outcomes,"black"),
       x=0.5, y=0.5, xjust=0.5, yjust=0.5, cex=2.5, bty="n")

dev.off()

###### 7. Plotting Figure S9 ######
#Preparing benefits data
aux.FigS9.BEN <- FigS9.BEN.Mi
aux.FigS9.BEN$aux.zero <- 0
aux.FigS9.BEN <- aux.FigS9.BEN[,c(3,1,2)]
aux.FigS9.BEN$SUM_NEI_BothOutcomes <- rowSums(aux.FigS9.BEN[,c(2,3)])
aux.FigS9.BEN$VMT <- DSB.EFs.NEI$VMT
aux.FigS9.BEN <- aux.FigS9.BEN[order(aux.FigS9.BEN$SUM_NEI_BothOutcomes,decreasing=TRUE),]
aux.FigS9.BEN$Cum.BEN <- aux.FigS9.BEN$SUM_NEI_BothOutcomes * aux.FigS9.BEN$VMT

aux.FigS9.BEN$VMT.Left <- 0
aux.FigS9.BEN$VMT.Right <- 0
aux.FigS9.BEN$VMT.Right[1] <- aux.FigS9.BEN$VMT[1]

aux.FigS9.BEN$Cum.BEN.Left <- 0
aux.FigS9.BEN$Cum.BEN.Right <- 0
aux.FigS9.BEN$Cum.BEN.Right[1] <- aux.FigS9.BEN$Cum.BEN[1]

for (i in 2:length(aux.FigS9.BEN$VMT))
{
  aux.FigS9.BEN$VMT.Left[i] <- sum(aux.FigS9.BEN$VMT[1:(i-1)])
  aux.FigS9.BEN$VMT.Right[i] <- sum(aux.FigS9.BEN$VMT[1:(i)])
  
  aux.FigS9.BEN$Cum.BEN.Left[i] <- sum(aux.FigS9.BEN$Cum.BEN[1:(i-1)])
  aux.FigS9.BEN$Cum.BEN.Right[i] <- sum(aux.FigS9.BEN$Cum.BEN[1:(i)])
}
Ben.Cols <- 1:3

#Preparing DSB impacts data
aux.FigS9.DSB <- FigS9.DSB.Mi
aux.FigS9.DSB$aux.zero <- 0
aux.FigS9.DSB <- aux.FigS9.DSB[,c(6,1,3,5,4,2)]

aux.FigS9.DSB$SUM_NEI_BothOutcomes <- rowSums(aux.FigS9.DSB[,c(2:6)])
aux.FigS9.DSB$VMT <- DSB.EFs.NEI$VMT
aux.FigS9.DSB <- aux.FigS9.DSB[order(aux.FigS9.DSB$SUM_NEI_BothOutcomes,decreasing=TRUE),]
aux.FigS9.DSB$Cum.DSB <- aux.FigS9.DSB$SUM_NEI_BothOutcomes * aux.FigS9.DSB$VMT

aux.FigS9.DSB$VMT.Left <- 0
aux.FigS9.DSB$VMT.Right <- 0
aux.FigS9.DSB$VMT.Right[1] <- aux.FigS9.DSB$VMT[1]

aux.FigS9.DSB$Cum.DSB.Left <- 0
aux.FigS9.DSB$Cum.DSB.Right <- 0
aux.FigS9.DSB$Cum.DSB.Right[1] <- aux.FigS9.DSB$Cum.DSB[1]

for (i in 2:length(aux.FigS9.DSB$VMT))
{
  aux.FigS9.DSB$VMT.Left[i] <- sum(aux.FigS9.DSB$VMT[1:(i-1)])
  aux.FigS9.DSB$VMT.Right[i] <- sum(aux.FigS9.DSB$VMT[1:(i)])
  
  aux.FigS9.DSB$Cum.DSB.Left[i] <- sum(aux.FigS9.DSB$Cum.DSB[1:(i-1)])
  aux.FigS9.DSB$Cum.DSB.Right[i] <- sum(aux.FigS9.DSB$Cum.DSB[1:(i)])
}
Impact.Cols <- 1:6

#Plotting
aux.pdffilename <- paste0(FigPath,"Fig_","S9","_",FigVersion,".pdf")
aux.right.axis.scaling <- (5/2)*1e-9

aux.rec.border <- c(rep(NA,5),"black")


pdf(aux.pdffilename,height=1.5*8.4,width=1.5*6.5)

aux.col.cum <- "red"
par(cex.lab=1,cex.main=1.25,cex.axis=0.8, 
    mar=c(5,3.5,3,3.5), mgp=c(2,0.5,0),
    mfcol=c(2,1))
plot(NA,xlim=c(0,7),ylim=c(0,6),
     xaxt="n",yaxt="n",
     xlab=expression(bold("Cumulative School Bus Miles Travelled [billion miles]")),
     ylab=expression(bold("Health Benefits per mile [2022 USD per mile]")),
     main=expression(bold("(a) Benefits of Electric School Buses: Fleet in 2017 [cents per mile]")))

for(j in 1:2)
{
  rect(xleft=aux.FigS9.BEN$VMT.Left*1e-9,
       xright=aux.FigS9.BEN$VMT.Right*1e-9,
       ybottom=rowSums(as.matrix(aux.FigS9.BEN[,Ben.Cols[1:(j)]])),
       ytop=rowSums(as.matrix(aux.FigS9.BEN[,Ben.Cols[1:(j+1)]])),
       border=aux.rec.border[j], lwd=0.001, col=col.outcomes[j])
}

summary(rowSums(as.matrix(aux.FigS9.BEN[,Ben.Cols[1:(j)]])))


segments(x0=0,
         x1=aux.FigS9.BEN$VMT.Right[1]*1e-9,
         y0=0,
         y1=aux.right.axis.scaling*aux.FigS9.BEN$Cum.BEN.Right[1],
         col=aux.col.cum, lwd=1, lty=2)

segments(x0=aux.FigS9.BEN$VMT.Left*1e-9,
         x1=aux.FigS9.BEN$VMT.Right*1e-9,
         y0=aux.right.axis.scaling*aux.FigS9.BEN$Cum.BEN.Left,
         y1=aux.right.axis.scaling*aux.FigS9.BEN$Cum.BEN.Right,
         col=aux.col.cum, lwd=2, lty=1)

aux.US.Mean <- sum(aux.FigS9.BEN$Cum.BEN)/sum(aux.FigS9.BEN$VMT)
segments(x0=-10,
         x1=sum(aux.FigS9.BEN$VMT)*1e-9,
         y0=aux.US.Mean,
         y1=aux.US.Mean,
         col="darkblue",lty=1,lwd=2)

par(mgp=c(2,0.5,0))
axis(side=2,at=c(0.5*0:12),
     labels=paste0("$",format(0.5*0:12,nsmall=2)),
     font=2,font.axis=2)


par(mgp=c(2,0.5,0))
axis(side=1,at=c(1*0:6,sum(aux.FigS9.BEN$VMT)*1e-9),labels=c(1*0:6,6.8),font=2,font.axis=2)

aux.legend.text <- c(expression(paste("Mortality")),
                     expression(paste("Childhood Asthma")),
                     expression(paste("U.S. VMT-weighted mean ($0.28/mile)")),
                     expression(paste("Cumulative Health Benefit")))

legend(x=0.5,y=6.2,
       xjust=0, yjust=1,
       legend=aux.legend.text,
       pch=c(rep(22,2),NA,NA),pt.cex=c(rep(1.2,2),NA,NA),col=c(rep(NA,2),"darkblue",aux.col.cum),
       pt.bg=c(col.outcomes[1:2],NA,NA),lty=c(rep(NA,2),1,1),lwd=c(rep(NA,2),2,2),
       box.lwd=NA, bg="white", cex=1)
box(which="plot",lty="solid")

par(mgp=c(2,0.5,0))

axis(side=4,at=c(0.5*0:12),
     labels=c("0","","$400M","","$800M","",
              "$1.2B","","$1.6B","","$2.0B","","$2.4B"),
     col.ticks = aux.col.cum, col.axis=aux.col.cum, col=aux.col.cum,font=2,font.axis=2)
mtext("Cumulative Health Benefit Per Year [2022 USD]",
      side=4,line=2,at=3, col=aux.col.cum, font=2)

par(new=T)

aux.col.cum <- "steelblue3"

par(cex.lab=1,cex.main=1.25,cex.axis=0.8, 
    mar=c(5,3.5,3,3.5), mgp=c(2,0.5,0),
    mfcol=c(2,1))
plot(NA,xlim=c(0,7),ylim=c(0,6),
     xaxt="n",yaxt="n",
     xlab=expression(bold("Cumulative School Bus Miles Travelled [billion miles]")),
     ylab=expression(bold("Health Impacts per mile [2022 USD per mile]")),
     main=expression(bold("(b) Diesel School Bus Health Impacts: 2017 Fleet [cents per mile]")))

for(j in 1:5)
{
  rect(xleft=aux.FigS9.DSB$VMT.Left*1e-9,
       xright=aux.FigS9.DSB$VMT.Right*1e-9,
       ybottom=rowSums(as.matrix(aux.FigS9.DSB[,Impact.Cols[1:(j)]])),
       ytop=rowSums(as.matrix(aux.FigS9.DSB[,Impact.Cols[1:(j+1)]])),
       border=aux.rec.border[j], lwd=0.001, col=col.pols[j])
}

summary(rowSums(as.matrix(aux.FigS9.DSB[,Impact.Cols[1:(j)]])))


segments(x0=0,
         x1=aux.FigS9.DSB$VMT.Right[1]*1e-9,
         y0=0,
         y1=aux.right.axis.scaling*aux.FigS9.DSB$Cum.DSB.Right[1],
         col=aux.col.cum, lwd=1, lty=2)

segments(x0=aux.FigS9.DSB$VMT.Left*1e-9,
         x1=aux.FigS9.DSB$VMT.Right*1e-9,
         y0=aux.right.axis.scaling*aux.FigS9.DSB$Cum.DSB.Left,
         y1=aux.right.axis.scaling*aux.FigS9.DSB$Cum.DSB.Right,
         col=aux.col.cum, lwd=2, lty=1)

aux.US.Mean <- sum(aux.FigS9.DSB$Cum.DSB)/sum(aux.FigS9.DSB$VMT)
segments(x0=-10,
         x1=sum(aux.FigS9.DSB$VMT)*1e-9,
         y0=aux.US.Mean,
         y1=aux.US.Mean,
         col="cyan",lty=1,lwd=2)

par(mgp=c(2,0.5,0))
axis(side=2,at=c(0.5*0:12),
     labels=paste0("$",format(0.5*0:12,nsmall=2)),
     font=2,font.axis=2)


par(mgp=c(2,0.5,0))
axis(side=1,at=c(1*0:6,sum(aux.FigS9.DSB$VMT)*1e-9),labels=c(1*0:6,6.8),font=2,font.axis=2)

aux.legend.text <- c(expression(paste("Primary PM"[2.5],"",sep="")),
                     expression(paste("SO"[2],"",sep="")),
                     expression(paste("NO"["x"],"",sep="")),
                     expression(paste("NH"[3],"",sep="")),
                     expression(paste("VOC", "",sep="")),
                     expression(paste("U.S. VMT-weighted mean ($0.29/mile)")),
                     expression(paste("Cumulative Health Impact")))
aux.pol.order <- c(1,3,5,4,2)

legend(x=0.5,y=6.2,
       xjust=0, yjust=1,
       legend=aux.legend.text[c(aux.pol.order,6,7)],
       pch=c(rep(22,5),NA,NA),pt.cex=c(rep(1.2,5),NA,NA),col=c(rep(NA,5),"cyan",aux.col.cum),
       pt.bg=c(aux.col.pol[aux.pol.order],NA,NA),lty=c(rep(NA,5),1,1),lwd=c(rep(NA,5),2,2),
       box.lwd=NA, bg="white", cex=1)
box(which="plot",lty="solid")

par(mgp=c(2,0.5,0))

axis(side=4,at=c(0.5*0:12),
     labels=c("0","","$400M","","$800M","",
              "$1.2B","","$1.6B","","$2.0B","","$2.4B"),
     col.ticks = aux.col.cum, col.axis=aux.col.cum, col=aux.col.cum,font=2,font.axis=2)
mtext("Cumulative Health Impact Per Year [2022 USD]",
      side=4,line=2,at=3, col=aux.col.cum, font=2)

dev.off()

###### 8. Mapping counties to school districts and creating supplementary datasets ######
#What it will contain: 
#SD.01 -- Results by COUNTY per MILE
#Benefits (BEN) -- by EF and outcome
#DSB health impacts (DSB) -- by EF, outcome, and pollutant

#SD.02 -- Results by COUNTY per BUS
#Benefits (BEN) -- by EF and outcome
#DSB health impacts (DSB) -- by EF, outcome, and pollutant

#SD.03 -- Results by SCHOOL DISTRICT per MILE
#Benefits (BEN) -- by EF and outcome
#DSB health impacts (DSB) -- by EF, outcome, and pollutant

#SD.04 -- Results by SCHOOL DISTRICT per BUS
#Benefits (BEN) -- by EF and outcome
#DSB health impacts (DSB) -- by EF, outcome, and pollutant

#SD.05 -- Marginal damages of ground-level emissions by COUNTY
#By pollutant and outcome

#SD.06 -- ESB impacts per mile and per bus, by pollutant

###### 8.0. Creating basic SD formats ######
#Creating basic SD format for counties
SD.County.Format <- as.data.frame(c(rep("County",3108),c("All Large Central Metro Counties","All Large Fringe Metro Counties",
                                              "All Medium Metro Counties","All Small Metro Counties",
                                              "All Micropolitan Counties","All Non-core ('Rural') Counties",
                                              "All Large Metro Counties (Central + Fringe)",
                                              "U.S. -- All Counties in the Contiguous United States")))
names(SD.County.Format) <- "Geography"
SD.County.Format$FIPS_STCOU <- c(DSB.EFs.NEI$FIPS_STCOU,rep(NA,8))
SD.County.Format$Row_Number <- 1:3116

aux.SD.County.Format <- subset(DSB.EFs.NEI, select=c("FIPS_STCOU","COUNTY_NAME","FIPS_STATE","STATE_NAME"))
SD.County.Format <- merge(x=SD.County.Format,y=aux.SD.County.Format,by="FIPS_STCOU",all=TRUE)

SD.County.Format <- merge(x=SD.County.Format,y=Urb,by.x="FIPS_STCOU",by.y="STCOU",all=TRUE)

names(SD.County.Format)
SD.County.Format <- subset(SD.County.Format, select=c("Row_Number","Geography",
                                        "FIPS_STCOU","COUNTY_NAME","FIPS_STATE","STATE_NAME","NCHS.Urb"))
  
names(SD.County.Format)[4:7] <- c("County_Name","FIPS_State","State_Name","NCHS_Urban_Rural_Classification")
SD.County.Format <- SD.County.Format[order(SD.County.Format$Row_Number,decreasing=FALSE),]

sum(SD.County.Format$FIPS_STCOU[1:3108]==DSB.EFs.NEI$FIPS_STCOU[1:3108])
sum(SD.County.Format$FIPS_STCOU[1:3108]!=DSB.EFs.NEI$FIPS_STCOU[1:3108])

#Creating basic school district
load("Inputs/SchoolDistrict_County_Mapping.RData")
School.Districts <- SchoolDistrict.County.Mapping
names(School.Districts)[6:7] <- c("FIPS_STCOU","Share_of_SchoolDistrict_in_STCOU")

#Restricting to contiguous U.S.
School.Districts <- School.Districts[!School.Districts$ST %in% c(2,15,72),] #AK, HI, PR

School.Districts[!School.Districts$FIPS_STCOU %in% DSB.EFs.NEI$FIPS_STCOU,]

SD.District.Format <- unique(School.Districts[,1:5])
length(unique(SD.District.Format$LEAID))

SD.District.Format <- SD.District.Format[order(SD.District.Format$LEAID,decreasing=FALSE),]
SD.District.Format$Row_Number <- 1:length(SD.District.Format$LEAID) 
names(SD.District.Format)
SD.District.Format <- subset(SD.District.Format, select=c("Row_Number","LEAID","NAME_LEA23","ST","STATE_NAME"))

###### 8.1. Supplemental Dataset 1 -- SD.01 ######
SD.01 <- SD.County.Format

aux.ef <- c("NEI","MY2005","MY2010","MY2020")
aux.oc <- c("BothOutcomes","Deaths","Asthma")


for (i in 1:length(aux.ef))
{
  aux.subset <- paste("SUM",aux.ef[i],aux.oc,sep="_")
  aux.subset2 <- subset(DSB.Mi, select=aux.subset)
  aux.subset3 <- subset(DSB.Mi.URB, select=aux.subset)
  aux.subset2 <- rbind(aux.subset2,aux.subset3)
  
  aux.b <- grep("BothOutcomes",colnames(aux.subset2))
  aux.subset2[,aux.b] <- aux.subset2[,aux.b] - ESB.Mi$SUM_BothOutcomes
  SD.01$newvar <- aux.subset2[,aux.b]
  names(SD.01)[dim(SD.01)[2]] <- paste(aux.ef[i],"BEN","Mile","SUM","AllOutcomes",sep="_")
  
  aux.d <- grep("Deaths",colnames(aux.subset2))
  aux.subset2[,aux.d] <- aux.subset2[,aux.d] - ESB.Mi$SUM_Deaths
  SD.01$newvar <- aux.subset2[,aux.d]
  names(SD.01)[dim(SD.01)[2]] <- paste(aux.ef[i],"BEN","Mile","SUM","Deaths",sep="_")
  
  aux.a <- grep("Asthma",colnames(aux.subset2))
  aux.subset2[,aux.a] <- aux.subset2[,aux.a] - ESB.Mi$SUM_Asthma
  SD.01$newvar <- aux.subset2[,aux.a]
  names(SD.01)[dim(SD.01)[2]] <- paste(aux.ef[i],"BEN","Mile","SUM","Asthma",sep="_")
}

aux.ef <- c("NEI","MY2005","MY2010","MY2020")
aux.pol.nei <- c("SUM","PM25","SO2","NOX","NH3","VOC")
aux.pol.greet <- aux.pol.nei[c(1,2,4,6)]
aux.oc <- c("BothOutcomes","Deaths","Asthma")
aux.oc.name <- c("AllOutcomes","Deaths","Asthma")

for (i in 1:length(aux.ef))
{
 aux.sel.pol <- aux.pol.greet
 if(i==1)
 {
   aux.sel.pol <- aux.pol.nei
 }
  for(j in 1:length(aux.sel.pol))
  {
    for(k in 1:length(aux.oc))
    {
      aux.subset <- paste(aux.sel.pol[j],aux.ef[i],aux.oc[k],sep="_")
      aux.filt <- which(names(DSB.Mi) == aux.subset)
      aux.subset2 <- as.numeric(DSB.Mi[,aux.filt])
      
      aux.filt <- which(names(DSB.Mi.URB) == aux.subset)
      aux.subset3 <- as.numeric(DSB.Mi.URB[,aux.filt])
      
      aux.subset2 <- c(aux.subset2,aux.subset3)

      print(paste(i,j,k,names(aux.subset2),dim(aux.subset)[2],sep="---"))
      SD.01$newvar <- aux.subset2
      names(SD.01)[dim(SD.01)[2]] <- paste(aux.ef[i],"DSB","Mile",aux.sel.pol[j],aux.oc.name[k],sep="_")
    }
  }
}

SD.01 <- SD.01[order(SD.01$Row_Number,decreasing=FALSE),]


###### 8.2. Supplemental Dataset 2 -- SD.02 ######
SD.02 <- SD.County.Format

aux.ef <- c("NEI","MY2005","MY2010","MY2020")
aux.oc <- c("BothOutcomes","Deaths","Asthma")


for (i in 1:length(aux.ef))
{
  aux.subset <- paste("SUM",aux.ef[i],aux.oc,sep="_")
  aux.subset2 <- subset(DSB.BUS, select=aux.subset)
  aux.subset3 <- subset(DSB.BUS.URB, select=aux.subset)
  aux.subset2 <- rbind(aux.subset2,aux.subset3)
  
  aux.b <- grep("BothOutcomes",colnames(aux.subset2))
  aux.subset2[,aux.b] <- aux.subset2[,aux.b] - ESB.BUS$SUM_BothOutcomes
  SD.02$newvar <- aux.subset2[,aux.b]
  names(SD.02)[dim(SD.02)[2]] <- paste(aux.ef[i],"BEN","BUS","SUM","AllOutcomes",sep="_")
  
  aux.d <- grep("Deaths",colnames(aux.subset2))
  aux.subset2[,aux.d] <- aux.subset2[,aux.d] - ESB.BUS$SUM_Deaths
  SD.02$newvar <- aux.subset2[,aux.d]
  names(SD.02)[dim(SD.02)[2]] <- paste(aux.ef[i],"BEN","BUS","SUM","Deaths",sep="_")
  
  aux.a <- grep("Asthma",colnames(aux.subset2))
  aux.subset2[,aux.a] <- aux.subset2[,aux.a] - ESB.BUS$SUM_Asthma
  SD.02$newvar <- aux.subset2[,aux.a]
  names(SD.02)[dim(SD.02)[2]] <- paste(aux.ef[i],"BEN","BUS","SUM","Asthma",sep="_")
}

aux.ef <- c("NEI","MY2005","MY2010","MY2020")
aux.pol.nei <- c("SUM","PM25","SO2","NOX","NH3","VOC")
aux.pol.greet <- aux.pol.nei[c(1,2,4,6)]
aux.oc <- c("BothOutcomes","Deaths","Asthma")
aux.oc.name <- c("AllOutcomes","Deaths","Asthma")


for (i in 1:length(aux.ef))
{
  aux.sel.pol <- aux.pol.greet
  if(i==1)
  {
    aux.sel.pol <- aux.pol.nei
  }
  for(j in 1:length(aux.sel.pol))
  {
    for(k in 1:length(aux.oc))
    {
      aux.subset <- paste(aux.sel.pol[j],aux.ef[i],aux.oc[k],sep="_")
      aux.filt <- which(names(DSB.BUS) == aux.subset)
      aux.subset2 <- as.numeric(DSB.BUS[,aux.filt])
      
      aux.filt <- which(names(DSB.BUS.URB) == aux.subset)
      aux.subset3 <- as.numeric(DSB.BUS.URB[,aux.filt])
      
      aux.subset2 <- c(aux.subset2,aux.subset3)

      print(paste(i,j,k,names(aux.subset2),dim(aux.subset)[2],sep="---"))
      SD.02$newvar <- aux.subset2
      names(SD.02)[dim(SD.02)[2]] <- paste(aux.ef[i],"DSB","BUS",aux.sel.pol[j],aux.oc.name[k],sep="_")
    }
  }
}

SD.02 <- SD.02[order(SD.02$Row_Number,decreasing=FALSE),]

###### 8.3. Supplemental Dataset 3 -- SD.03 ######
#Creating basic SD format for school districts
SD.03 <- SD.District.Format

aux.variables <- names(SD.01)[8:73]
aux.sd01 <- subset(SD.01,select=c("FIPS_STCOU",aux.variables))
aux.sd01 <- aux.sd01[!is.na(aux.sd01$FIPS_STCOU),]
i<-1
for(i in 1:length(aux.variables))
{
  aux.var1 <- aux.variables[i]
  aux.sel1 <- subset(aux.sd01,select=c("FIPS_STCOU",aux.var1))
  aux.dis1 <- merge(x=School.Districts,y=aux.sel1,by="FIPS_STCOU",all.x=TRUE)
  aux.which1 <- which(names(aux.dis1)==aux.var1)
  aux.dis1$newvar <- aux.dis1$Share_of_SchoolDistrict_in_STCOU * aux.dis1[,aux.which1]
  aux.agg1 <- aggregate(aux.dis1$newvar, by=list(aux.dis1$LEAID), FUN=sum)
  aux.test <- sum(names(aux.agg1) == c("Group.1","x"))
  if(aux.test != 2)
  {
    print(paste("ERROR",i,sep="---"))
  }
  names(aux.agg1)[names(aux.agg1)=="Group.1"] <- "LEAID"
  names(aux.agg1)[names(aux.agg1)=="x"] <- "newvar"
  SD.03 <- merge(x=SD.03,y=aux.agg1,by="LEAID",all=TRUE)
  aux.which3 <- which(names(SD.03)=="newvar")
  names(SD.03)[aux.which3] <- aux.var1
}

sum(is.na(SD.03))
SD.03 <- SD.03[order(SD.03$Row_Number,decreasing=FALSE),]
names(SD.03)
SD.03 <- subset(SD.03,select=c(names(SD.03)[c(2,1,3:dim(SD.03)[2])]))

School.Districts[grep("New York",School.Districts$NAME_LEA23),]


###### 8.4. Supplemental Dataset 4 -- SD.04 ######
#Creating basic SD format for school districts
SD.04 <- SD.District.Format

aux.variables <- names(SD.02)[8:73]
aux.sd02 <- subset(SD.02,select=c("FIPS_STCOU",aux.variables))
aux.sd02 <- aux.sd02[!is.na(aux.sd02$FIPS_STCOU),]

for(i in 1:length(aux.variables))
{
  aux.var1 <- aux.variables[i]
  aux.sel1 <- subset(aux.sd02,select=c("FIPS_STCOU",aux.var1))
  aux.dis1 <- merge(x=School.Districts,y=aux.sel1,by="FIPS_STCOU",all.x=TRUE)
  aux.which1 <- which(names(aux.dis1)==aux.var1)
  aux.dis1$newvar <- aux.dis1$Share_of_SchoolDistrict_in_STCOU * aux.dis1[,aux.which1]
  aux.agg1 <- aggregate(aux.dis1$newvar, by=list(aux.dis1$LEAID), FUN=sum)
  aux.test <- sum(names(aux.agg1) == c("Group.1","x"))
  if(aux.test != 2)
  {
    print(paste("ERROR",i,sep="---"))
  }
  names(aux.agg1)[names(aux.agg1)=="Group.1"] <- "LEAID"
  names(aux.agg1)[names(aux.agg1)=="x"] <- "newvar"
  SD.04 <- merge(x=SD.04,y=aux.agg1,by="LEAID",all=TRUE)
  aux.which3 <- which(names(SD.04)=="newvar")
  names(SD.04)[aux.which3] <- aux.var1
}

sum(is.na(SD.04))

SD.04 <- SD.04[order(SD.04$Row_Number,decreasing=FALSE),]
names(SD.04)
SD.04 <- subset(SD.04,select=c(names(SD.04)[c(2,1,3:dim(SD.04)[2])]))

summary(as.numeric(as.matrix(SD.04[,-c(1:5)]))/as.numeric(as.matrix(SD.03[,-c(1:5)])))

School.Districts[grep("New York",School.Districts$NAME_LEA23),]

###### 8.5. Supplemental Dataset 5 -- SD.05 ######
SD.05 <- SD.County.Format
SD.05 <- SD.05[!is.na(SD.05$FIPS_STCOU),]


sum(order(SD.05$FIPS_STCOU,decreasing=FALSE)==1:3108)
sum(order(SD.05$FIPS_STCOU,decreasing=FALSE)!=1:3108)

sum(SD.05$FIPS_STCOU==SD.05$FIPS_STCOU[order(SD.05$FIPS_STCOU,decreasing=FALSE)])
sum(SD.05$FIPS_STCOU!=SD.05$FIPS_STCOU[order(SD.05$FIPS_STCOU,decreasing=FALSE)])


aux.sd05 <- cbind(Marg.Damages.GEMM_BasePM_2017_Mortality_2017,Marg.Damages.Tetreault_BasePM_2017_Asthma_2019[,c(16:20)])
#Sum
aux.sum <- aux.sd05[,1:5] + aux.sd05[,6:10]

aux.sd05 <- cbind(aux.sum,aux.sd05)
aux.eg <- expand.grid(c("PM25","SO2","NOX","NH3","VOC"),c("BothOutcomes","Deaths","Asthma"))
colnames(aux.sd05) <- paste("MargDamages",aux.eg[,1],aux.eg[,2],sep="_")

SD.05 <- cbind(SD.05,aux.sd05)

SD.05 <- SD.05[order(SD.05$Row_Number,decreasing=FALSE),]


###### 8.6. Supplemental Dataset 6 -- SD.06 ######
names(ESB.Mi)
ESB.Mi

SD.01[3116,]

aux.ESB <- ESB.Mi[,-(grep("Attrib",names(ESB.Mi)))]
aux.ESB.bus <- aux.ESB * School.Bus.LT.Miles
names(aux.ESB) <- paste0("ESB_Mile_",names(aux.ESB))
names(aux.ESB.bus) <- paste0("ESB_BUS_",names(aux.ESB.bus))
aux.ESB.3 <- aux.ESB
aux.ESB.3$Geography <- "U.S. Average"
aux.ESB.3 <- subset(aux.ESB.3, select="Geography")

SD.06 <- cbind(aux.ESB.3, aux.ESB, aux.ESB.bus)


write.csv(SD.01,file="Outputs/Supplemental_Datasets/SD01.csv",row.names=FALSE)
write.csv(SD.02,file="Outputs/Supplemental_Datasets/SD02.csv",row.names=FALSE)
write.csv(SD.03,file="Outputs/Supplemental_Datasets/SD03.csv",row.names=FALSE)
write.csv(SD.04,file="Outputs/Supplemental_Datasets/SD04.csv",row.names=FALSE)
write.csv(SD.05,file="Outputs/Supplemental_Datasets/SD05.csv",row.names=FALSE)
write.csv(SD.06,file="Outputs/Supplemental_Datasets/SD06.csv",row.names=FALSE)

###### 9. Plotting Fig 2 ######
library("sp")
library("raster")
library("viridis")

#Loading U.S. Census Shapefile with county boundaries to plot the map in Figure 2
# Source: U.S. Census Bureau, 2019 Cartographic Boundaries: County: 500k resolution level. 
# Available at: https://www2.census.gov/geo/tiger/GENZ2019/shp/cb_2019_us_county_500k.zip (Accessed 21 September 2019).
Counties <- shapefile("Inputs/USCB_County_Shapefile_500k/cb_2019_us_county_500k.shp")

Counties$STCOU <- as.numeric(as.character(Counties$GEOID))
Counties.ContUS <- Counties[Counties$STCOU %in% DSB.EFs.NEI$FIPS_STCOU,]

aux.data <- subset(DSB.EFs.NEI,select=c("FIPS_STCOU"))
aux.data$BEN.BUS <- Fig2.BEN.BUS

Counties.ContUS <- merge(x=Counties.ContUS,y=aux.data, 
                         by.x="STCOU", by.y="FIPS_STCOU",
                         all.x=TRUE, all.y=TRUE)

Counties.ContUS[Counties.ContUS$BEN.BUS<0] <- 1 #
Counties.ContUS$log10.BEN.BUS <- log10(Counties.ContUS$BEN.BUS)
Counties.ContUS$log10.BEN.BUS[Counties.ContUS$log10.BEN.BUS<=4] <- 4 #Benefits <= 10,000
summary(Counties.ContUS$log10.BEN.BUS)

#Creating color scale
aux.breaks <- c(3.99+0.01*(0:193)) 
col.scale <- viridis(n=(length(aux.breaks)+1),alpha=1,begin=0,end=1,direction=1)

aux.leg.at <- c(log10(c(10000,20000,30000,50000,100000,200000,300000,500000)),max(Counties.ContUS$log10.BEN.BUS))
#What the labels will indicate on the ticks
aux.leg.labels <- c("<10k","20k","30k","50k","100k","200k","300k","500k","820k")

aux.width <- 9*1.4
aux.height <- 5*1.4

plot(0)

aux.pdffilename <- paste0(FigPath,"Fig_","2","_",FigVersion,".pdf")

pdf(aux.pdffilename,
    width=aux.width,
    height=aux.height)

par(cex.lab=1.5)

spplot(Counties.ContUS, c("log10.BEN.BUS"), lwd=NA,
       col.regions = col.scale, at=aux.breaks,
       ylim=c(Counties.ContUS@bbox[2,1],(Counties.ContUS@bbox[2,2])),
       xlim=Counties.ContUS@bbox[1,],
       strip=FALSE,
       
       main=list(label=expression(bold(paste("Benefits of electric buses -- replacing 2017 Fleet Average"))), cex=1.5),
       xlab=list(label=expression(bold("Benefits of electric bus adoption in each county (2022 USD per bus)")), cex=1.5),
       colorkey=list(height=1, width=3, labels=list(labels=aux.leg.labels,at=aux.leg.at,cex=1,font.axis=2)))
dev.off()